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Abstract. We present a review of the quantum three-body problem, with emphasis 
on the different methodologies, different three-body atomic systems and their historical 
interest. With the review as the background, a more recently proposed non- variational, 
kinetic energy operator approach to the solution of quantum three-body problem 
is presented, based on the utilization of symmetries intrinsic to the kinetic energy 
operator, i.e., the three-body Laplacian operator with the respective masses. Through 
a four-step reduction process, the nine dimensional problem is reduced to a one 
dimensional coupled system of ordinary differential equations, amenable to accurate 
numerical solution as an infinite-dimensional algebraic eigenvalue problem. A key 
observation in this reduction process is that in the functional subspacc of the kinetic 
energy operator where all the rotational degrees of freedom have been projected out, 
there is an intrinsic symmetry which can be made explicit through the introduction of 
Jacobi-spherical coordinates. A numerical scheme is presented whereby the Coulomb 
matrix elements are calculated to a high degree of accuracy with minimal effort, and 
the truncation of the linear equations is carried out through a systematic procedure. 
The resulting matrix equations are solved through an iteration process. Numerical 
results are presented for (1) the negative hydrogen ion H _ , (2) the helium and helium- 
like ions [Z = 3 ~ 6), (3) the hydrogen molecular ion , and (4) the positronium 
negative ion Ps _ . Up to thirteen-significant-figure accuracy is achieved for the ground 
state eigenvalues when double precision programming is used. Comparison with the 
variational and other approaches shows our ground state eigenvalues to be comparable, 
generally with less decimal digits than the variational results, but can yield highly 
accurate wavefunctions as by-products. Results on low-lying excited states and their 
wavefunctions are obtained simultaneously with the ground state properties, some 
at accuracies not achieved by other methods. In particular, for the doubly excited 
state 3 P e of H~ and the 1 ^P e states of helium, some results are obtained for the 
first time. Also, we have calculated fourteen H^" excited states, up to its dissociation 
level. Analysis of the wavefunction characteristics, especially in relation to the electron- 
electron correlation effects, are presented. A significant advantage of the kinetic energy 
operator approach is its general applicability to different three-body systems, with only 
the charges, masses, and the symmetry of the desired state as the required inputs. 
Potential applications of the present approach to scattering and other problems are 
noted. 
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1. Introduction 

The dynamics of three interacting bodies constitutes one of the oldest challenges in 
physics. Studies in this field can be traced back to the work of Euler in the 18th century. 
In the beginning of the 20th century, the failure of the Bohr-Sommerfeld quantization 
to correctly describe the ground state of helium has stimulated the development of the 
"new" quantum theory, formulated by Schrodinger and Heisenberg. Almost 100 years 
after the founding of modern quantum mechanics, there is still a continuing effort to 
improve the solution methods or to invent new approaches for the seemingly simple 
three-body Coulomb system. 

In the early calculations of the two-electron systems, the focus was usually on the 
bound-state spectra, with the helium and helium-like ions as the proving ground. The 
spectra could be calculated efficiently with the help of the Hartree-Fock self-consistent- 
field method based on the variational principle. Very high accuracy can be achieved 
for eigenvalues such that they may be compared with high precision measurements. 
Bether and Salpeter P have summarized the early works in this area. The more recent 
calculations of the bound states have extended such routines with new numerical schemes 
and judicious choices of basis functions. 

In a seminal experiment by Madden and Codding j2j, the discovery of strong 
electron-electron correlation effects in doubly excited resonance states of helium has 
triggered the development of group-theoretical and adiabatic quantum approximations 
to understand these effects. The two-electron dynamics were again at the forefront of a 
revival. As doubly excited resonant states could not be tracked by an effective single- 
particle method, over the past four decades the effort to understand doubly excited 
resonances has stimulated much of the theoretical research on two-electron atoms. The 
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role of electron correlation has become more important with increasing order of double 
excitations, and the correlated three-body Coulomb dynamics have been found to cause 
an extremely rich and complicated resonance spectrum. Hence two-electron atoms have 
come to represent a prototype of the few-body systems strongly affected by electronic 
correlation. 

A more difficult problem is the three-body scattering, such as hydrogen-electron 
scattering. This problem has attracted much attention recently. By using the finite 
element method, Levin and Shertzer [Hj analyzed the S-wave phase shifts for low- 
energy positron-hydrogen scattering. Botero and Shertzer [3] directly solved the 
Schrodinger equation for electron-hydrogen scattering, and Rescigno et al j^j and 
Baertschy et al |B] used supercomputers to obtain a complete numerical solution of 
the hydrogen atom ionization through electron collision. 

Three-body systems, especially the two-electron systems, remain an active field of 
research today. This persistent interest can be traced to the fact that the three-body 
problem is just complex enough for rather sophisticated theoretical concepts, yet simple 
enough to provide accurate numerical and experimental tests. 

More recently, Hsiang and Hsiang [3 IH1 El QHj have outlined a new approach to the 
quantum three-body problem which was based on the systematic exploitation of all the 
intrinsic dynamic symmetries of the three-body kinetic energy operator (the three-body 
Laplacian with the respective masses), some of which not fully recognized previously. 
The purposes of this work are to implement this new mathematical formulation and 
to compare the present approach with the conventional variational approach in terms 
of numerical results for both eigenvalues and eigenfunctions of a number of quantum 
three-body systems. 

The main conclusions of this work are that the present approach offers not only 
systematic computability, requiring minimal ad hoc inputs in the computational process, 
but also achieves numerical accuracy for both the eigenvalues and eigenfunctions. The 
latter advantage is particularly significant for excited states. In particular, the wave 
function characteristics of excited states with strong electron-electron correlation can 
be accurately delineated. 

In what follows, we first review the various approaches to the quantum three-body 
problem in section followed by a detailed presentation of the kinetic energy operator 
approach in section El The formulation essentially consists of a four-step reduction 
process, in which the dynamic symmetries intrinsic to the kinetic energy operator are 
fully utilized. The end result of the reduction is a one-dimensional coupled system of 
ordinary differential equations, solvable as an infinite-dimensional algebraic eigenvalue 
problem. In section HJ a scheme is presented for the numerical implementation of our 
approach. Due to the high precision required for the Coulomb matrix elements, a special 
integration technique is used to evaluate both the matrix elements as well as the product 
of the potential energy matrix with a vector. A sparse matrix solver is then applied to 
solve the linear equations iteratively. Truncation of the infinite linear system is carried 
out by following a rule implied by the asymptotic behavior of the eigenvalues, leading 
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to sequences of numerical data from which one can apply a systematic extrapolation 
procedure. The results for some typical Coulomb systems are presented in section 
with comparisons to those obtained via other approaches. We also explore the properties 
of three-body wave functions, and discuss some of their physical significances. It should 
be noted that all the results were obtained by using the same program. Inputs are the 
symmetry of the state (i.e., S, P, D, or F), the mass ratios, and the sign of the charges. 
The article concludes by noting some potential applications of the present approach. 



2. Review of the various approaches 

2.1. Variational method 

Quantum variational method is most suitable for obtaining accurate results for the 
ground state or low-lying states, and in this regard it is superior to the perturbation 
methods. The basic idea of the variational method, sometime also denoted as the Ritz 
variational method, named after the pioneer of the approach, is to write the trial wave 
function ty tr in some arbitrarily chosen mathematical form with variational parameters, 

% r = % r (a,/3,-f,---), (1) 

and then adjust the parameters to obtain the minimum energy 

fV* tr H* tr dT 

through the solution of a system of coupled equations: 








dE tr \ 


[<*,P,7, ■ ■ ■) 




da 


dE tr \ 


>,/3,7r--) 




dp 


dE tr \ 


[a, /3,7, • • ■) 





; (3) 

It should be noted that the trial wave function(s) must satisfy the symmetry condition 
imposed by the Fermi-Dirac statistics. The solution to (J3J) yields minima of energy in 
the multidimensional parameter space. The lowest energy minimum is treated as the 
approximate ground state eigenvalue. As an example, a simple trial wave function for 
the ground state of helium-like atoms is ip(ri, r 2 , r 12 ) = exp(— (Z — o-)(ri +r 2 )), where a 
represents the screening effect in an approximate way. Minimizing the energy functional 
gives E = — (Z — 5/16) 2 , i.e., E=-2.85 a.u. for the helium ground state jT]. However, 
for H this trial wave function is noted to fail in obtaining a bounded ground state. 
For the excited states, the Hyllerass-Undheim theorem states that the remaining energy 
minima of Q, A 2 , A3, • • •, are also upper bounds to the exact eigenvalues E 2 , E 3 , ■ ■ •, 
provided that the spectrum is bounded from below. However, the calculation of precise 
excited state energies is more difficult compared to that of the ground state, due to 
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the appearance of the subsidiary condition that the eigenfunction of every excited state 
must be orthogonal to the eigenfunctions of all the lower-order states. This condition 
reduces considerably the number of available trial functions which may be chosen to 
approximate the eigenfunction (to be inserted in the variation integral) of the particular 
excited state. As a result, convergence is not nearly as good. 

Historically, Kellner ^1] was the first to use the variational principle to obtain a 
rather precise ground state energy, E=-2.895 a.u.. His results were improved upon by 
the variational calculations of Hylleraas ^21 CSj who obtained E=-2.9037 a.u. using 
a trial wave function with 38 variational parameters. This method was later used by 
Kinoshita fl] in large-scale variational calculations. In a bold move, Frankowski and 
Pekeris ^3] used more than 200 parameters in trial functions to obtain the energies 
of helium-like systems that were not surpassed for almost two decades. With the 
appearance of powerful modern computers, however, people can now include more 
than one thousand parameters in trial wave functions in order to obtain high accuracy 
results for the ground state of helium, hydrogen-like ions and some muonic molecular 

ions ma im mm Ha eh. 

As a modified version of the variational method, the complex rotation method, 
based on the dilatation analytic continuation |2"2"1 |2"3*1 123] . was extensively used to 
calculate the doubly excited states of the two-electron systems [213 1213 12H I2H I2H| • The 
basic idea is that after a complex rotation, r —>■ re 10 , is applied to the radial coordinate, 
the resulting Hamiltonian becomes complex, i.e., H{0) = e~ 2l9 T + e~ l6 V , where T is the 
kinetic energy operator and V the Coulomb potential. The complex Hamiltonian will 
yield complex eigenenergies as a result of applying the Ritz variation, in which the real 
part would correspond to the position of the doubly excited state, and the imaginary 
part its life-time. 

2.2. Hyperspherical coordinates method 

Hyperspherical coordinates were first introduced into atomic physics by Gronwall [HD] 
to study the analytic structure of the Schrodinger equation for the ground state of 
helium atom. The basic idea of the hyperspherical approach to the three-body systems 
is to express the two relative (to center of mass) coordinates as a single six-dimensional 
vector, and the nonrelativistic Schrodinger equation in the six-dimensional space is to 
be solved without reference to the "wave functions" associated with individual particles. 
The condition of particle exchange symmetry then becomes a boundary condition on 
the three-body wave function on the hypersurface. 

The hyperspherical approach has been applied to solve bound states and scattering 
problems in many different fields of physics and chemistry. Many of the earlier 
works dealt with the basic structure of the mathematical functions encountered in 
hyperspherical coordinates. Here we introduce the hyperspherical coordinates for two- 
electron atomic systems such as helium and hydrogen negative ions where the mass of 
the nucleus is treated as infinite. The hyperspherical coordinates are then obtained by 
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defining 

P = \J r i+ r l; a = arctan(r 2 /ri), (4) 

where p is the hyperradius which measures the size of the system, and a is the 
hyperangle. Thus the two vectors ri and r 2 are replaced by six coordinates (p,fl), 
where Q = (a, 9x, 0i, 9 2 , 2 ) denotes collectively the five angles, with (6i,if)i) being the 
spherical angles of electron i. In hyper spherical coordinates the two-electron equation 
is given by 

Ia 1 -Ia 2 ---- + — -e) f(n, r 2 ) = 0. (5) 



2 2 r 1 r 2 r 12 
After eliminating the first-order derivatives in the differential operators by expressing 

tf(ri,r 2 ) = i>{p,tt)/(p 5/2 sin a cos a), (6) 

an equation in terms of ijj(p, Q) is obtained: 

/ f)2 \2 \ 

- + - + 2E 4fi) = 0, (7) 



where 



dp 2 p 2 p 



a 2 ( d 2 \ 2 I 2 \ 1 

V oa z cos z a sin a / 4 



is the grand angular momentum operator, lj being the angular momentum operator for 
electron i and C/p is the total Coulomb interaction potential among the three charged 
particles, with C given by 

C(a, 6) = — + 1 ==. (9) 

cos a sin a yl — sin 2a cosy 12 

Here 812 is the angle between the two electrons with respect to the nucleus (as the 
origin). 

The straightforward solution approach is to expand ip(p, fl) by the eigenfunctions 
of A 2 , called hyperspherical harmonics. This method has been applied by a number of 
authors to H~ and He systems jSH E21 EM Oil, but the rate of convergence is rather 
slow. To improve the rate of convergence, Haftel and Mandelzweig [23 EH E3 OH OH] 
introduced an exponential factor in the expansion, if; = x$, where x is chosen to be 
of the form \ — exp{— a(ri + r 2 ) + 6r 12 }, with a and b to be obtained variationally or 
by some ansatz. If a and b are appropriately chosen, the singularity in the Coulomb 
potential can be removed, and it is then possible to expand $ in terms of hyperspherical 
harmonics with rapid convergence. Another common approach, the adiabatic expansion, 
was introduced by Fano and first applied by Macek [J^. Details of this method can be 
found in the review article by Fano [Hi] and a relevant book 
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2. 3. Perturbation, Hartree-Fock, and the finite element methods 

In the perturbation method, the Hamiltonian is split into two parts, H = Hq + XHi, 
where the perturbation term \H\ is small in a relative sense. For the Hamiltonian 
equation 

(H + XHi - E)V = 0, (10) 

the eigenfunction \l/ and eigenvalue E are expanded in powers of the small parameter 
A, namely 

oo oo 

E=Y / ^ n En, * = (11) 

n=0 n=0 

Substitution of these expansions into the Schrodinger equation and equating the 
coefficients for each power of A to zero lead to an infinite set of coupled linear equations: 

#o*o - £o*o = (12) 

#o*i + #i*o - £o*i - £i*o = (13) 



# *n + #l*n-l " E Em^n-m = 0. (14) 
m=0 

If we consider \l/o and Eq as known, we can obtains the first order perturbation energy 

E 1 = J %H^ dr. (15) 

For the calculation of E 2 , E 3 , ■ • •, more efforts are needed. There is no theorem to 
clearly tell one how to choose the perturbation A#i, and there can be different choices 
for the same problem. For large Z in helium-like ions, the interaction I/V12 between 
the electrons can be treated as the perturbation; for the excited state, there is the 
unsymmetrical choice [T], #0 = T + V\ + V2 and A#i = W, where 

V 1 {r 1 ) = --, V 2 (r 2 ) = - — , W = — --. (16) 
ri r 2 ri2 r 2 

In the Hartree-Fock method, the wave functions for the two-electron atom must 

obey overall antisymmetry. That means 

* = -^^i(r 1 )^(r 2 ) ± ^(r 2 )^(ri)], (17) 

where the + sign in the above equation is used when the spin state has antisymmetry. 
The form of the trial wave functions ^i(r) or ^2( r ) is unknown, but the aim is to find the 
most accurate form possible for the two functions ipi(r) or ^2(1") that would minimize the 
expectation value of the Hamiltonian, which is regarded as a functional of the two trial 
wave functions. The application of the general variational principle leads to two coupled 
differential-integral equations (the Euler-Lagrange equations). In essence, the Hartree- 
Fock method is a mean-field theory. It differs from the Ritz variational approach in 
that the variation in the Ritz method is carried out by using parameters, whereas in the 
Hartree-Fock method the variation is carried out by solving coupled integral-differential 
equations. 
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Recently, the finite element method (FEM) has been used to directly solve the 
Schrodinger equation, especially for the S state of some systems. The domain of the 
wave function is segmented into tetrahedrons, each serving as the domain of a local 
polynomial basis set. The approximation of the wave function as a linear combination 
of these local polynomials is called a finite-element description. FEM treatment for 
the helium atom in the infinite nuclear mass approximation HH1 021 ED] , and 

the hydrogen molecular ion in the Born-Oppenheimer approximation, have been 
presented by several authors jHH E21 EH1 EH E3 EE] • An adaptable FEM approach was 
used by Ackermann j^Z] and co-workers to obtain the energy values to a precision of 
10 -11 with moderate computational effort. The advantage of FEM lies in its flexibility 
because of the local basis, but the drawback is that one often has to face huge sparse 
matrices. 



3. Formulation 



In a three-body Coulomb system, the position vectors and the masses of the three 
particles are denoted by i\, and Mj, j — 1, 2, 3. The relative masses are defined as rrij = 
Mj/M, where M is the total mass, M = X)Mj, and J2 m j — 1- Schrodinger equation is 
given by 

-— A^> + = EV, 
2M 

3 A ■ 

A = £-> (18) 

where Aj is the Laplacian operator with respect to the j-th particle, A is defined as 
the kinetic energy operator, and V is the Coulomb potential of the three interacting 
particles: 

Z^Z? Z'iZ'i Zi Z 1 ) 

V = - — - J - 1 + ] — LJl t+i — LJ ^- (19) 

I r 2 — r 3| | r 3 — r l| l r l — r 2\ 

Here Zj is the electric charge of the j-th particle. In the above and in what follows, 
the length unit is the Bohr radius, h 2 /m e e 2 , m e is the electron mass, also the mass 
unit, and m e e 2 /h 2 is the energy unit. The above system is uniquely determined by 
the six parameters, {rrij, Zj, j = 1, 2, 3}. The Coulomb potential depends only on 
the distance between each pair of particles, and the system is invariant under spatial 
translation, rotation, and inversion. Furthermore, if the system consists of identical 
particles, invariance with respect to the permutation of identical particles must be 
imposed. 



3.1. Overview of the approach 

As implied by the name, the focus of the present approach is on the kinetic energy 
operator A. Before delving into the details, it would be helpful to give an overview 
of the formulation, which in essence consists of four reduction steps. In the first 
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step (section l3.2|) . we reduce the nine dimensional problem to a six dimensional problem 
by simply using coordinates relative to the center of mass. Jacobi vectors will be 
introduced at this stage to facilitate later developments. In the second step ( section l3~3"|) . 
angular momentum eigenf unctions will be presented which are in the null space of 
the kinetic energy operator. That is, when the angular momentum eigenf unctions 
are operated on by the kinetic energy operator, the result is zero. Thus when the 
total wave function is expanded in terms of the angular momentum eigenfunctions, 
the coefficients of the expansion, which are rotationally invariant, satisfy a reduced 
Schrodinger equation which is three-dimensional. The set of functions that satisfy this 
reduced Schrodinger equation constitutes a subspace in which all the rotational degrees 
of freedom have been projected out. In the third step of the reduction ( section I3.4J) . 
this 3D reduced Schrodinger equation is expressed in terms of the Jacobi-spherical 
coordinates, leading to a form with an angular operator whose eigenfunctions are the 
Jacobi polynomials (hence the denotation of Jacobi-spherical coordinate). It should be 
emphasized that the Jacobi-spherical symmetry is particular only to the kinetic energy 
operator and not to the Coulomb potential, in contrast to the rotational symmetry. 
Hence the total Hamiltonian does not have this symmetry. However, the Coulomb 
potential becomes separable in the Jacobi-spherical coordinates. That is, the Coulomb 
potential can be expressed as a product of two terms, one of which depends only on the 
radial coordinate of the Jacobi-spherical coordinates. The wave functions of this reduced 
Schrodinger equation can therefore be expanded in terms of the Jacobi polynomials with 
coefficients depending only on the radial coordinate of the Jacobi-spherical coordinates. 
In the fourth and last reduction step (section l3.5|) . the substitution of the this expansion 
into the reduced Schrodinger equation leads to a set of coupled ordinary differential 
equations (ODEs) which is now only one dimensional, i.e., the solution depends only 
on the radial coordinate. The solution to the set of ODEs can be expanded in terms of 
the Laguerre polynomials, and the coefficients of this expansion satisfy an infinite linear 
system which is amenable to numerical solution fsection l3.6J) . 

3.2. Jacobi vectors and coordinate transformations 

The configuration space of a given three-body system is a nine-dimensional space 
consisting of a triplet of position vectors, 



where rj is the position vector of j-th particle. Without the loss of generality, one can 
assume that the center of gravity is fixed at the origin, thus reducing the configuration 
space to the following six-dimensional reduced configuration space: 



3? 9 



{ri, r 2 , r 3 } 



(20) 




(21) 
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Following Jacobi, we introduce the kinematic metric on the reduced configuration space 
by defining the inner product as 

({aA,{b,}): >>,a ; .|, r (22) 
i=i 

This definition is related to the Lagrange's least-action principle, which was reformulated 
by Jacobi, leading to the geometric explanation of mechanics. The metric ds 2 is defined 
in terms of the kinetic energy T by setting 

ds 2 = ^dt 2 . (23) 
For example, in the case of an n-body system 

n 

ds 2 = Yl m i ( dx * + d y* + dz ") » ( 24 ) 

i=l 

where M is the total mass, and rrtj is the percentage of mass for the i-th 
particle. Lagrange's least-action principle of classical mechanics has a simple geometric 
explanation: trajectories of a given mechanical system are exactly those geodesic curves 
in the metric space of ds 2 , defined in the configuration space. 

To each given m-triangle, i.e., triplet {ij} with J2 m j r j — 0, we define a pair of 
Jacobi vectors x = {x\, x 2 , x 3 ) and y = (y\, y 2 , 2/3) in the reduced space 9ft 6 , given by 



/ m l / m 2 m 3 

x =i/i (m 2 (ri -r 2 ) +m 3 (ri -r 3 )), y=Jz (r 2 - r 3 ) . (25) 

y 1 — rrii V mi 

It should be noted here that in this work, indices 2, 3 are used to label identical particles 
in our three-body system (e.g., the two electrons in the two-electron systems). Hence 
the exchange antisymmetry of the Fermi-Dirac statistics is manifest in letting y — > — y. 
The reduced configuration space can be represented by 3ft 3 © 3ft 3 = {(x, y); x, y G 3ft 3 }, 
and the metric ds 2 , becomes 

ds 2 = dx • dx + dy • dy, (26) 
with the kinetic energy operator given by 

A = A x + A y . (27) 
Also, the total angular momentum operator is given by 

L = L X + Ly 

= - i (x x V x + y x V y ) . (28) 

By using the Jacobi vectors x and y, the motion of the center-of-mass is removed from 
()18jl . leading to a six-dimensional equation: 

-±{A X + A y )* + V9 = EV. (29) 
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3.2.1. Rotationally invariant variables For given two vectors x and y, we can construct 
three rotationally invariant polynomials, 

fx = x ■ x, f 2 = y ■ y, / 3 = x • y. (30) 

Any rotationally invariant quantity with respect to x and y should be a function of 
/2> Z^)- Since the metric ds 2 is rotationally invariant, thus it can be expressed as 

ds 2 = d/.,-, (31) 

where (gtj) is the inverse matrix of the following (g iJ ): 

{g ij ) = (V/i • V/i)- (32) 

The three variables (/i, / 2 , Z^) will be used to express the reduced Schrodinger equation 
below. 

3.2.2. Jacobi spherical coordinates The Hamiltonian of the three-body problem (Lapla- 
cian plus Coulomb potential) has space-rotation symmetry, but the symmetry of the 
kinetic energy operator is higher than that of the potential energy (Coulomb interac- 
tion). Moreover, there is a direct relation between the kinetic energy operator and 
the kinetic metric ds 2 . It is obvious that {fx,f2ifz) is not orthogonal; hence terms 
{g^ dfidfj,i 7^ j} do not vanish in ds 2 . By making use of the following coordinate 
transformation QD] , namely 



P = V/i + /a. cos ^ = 2(/i/a - /a) 1/2 (/i + h)~\ 
tan/3 = 2/ 3 (/ 2 - /1)" 1 , sin/3 = 2/ 3 [(/ 2 - /1) 2 + 4/ 2 



-1/2 



'3 

(0<p<cx), 0<^< tt/2, -tt < /3 < tt), (33) 
the metric ds 2 becomes 

ds 2 = dp 2 + ^p 2 (d# 2 + sin 2 6d{3 2 ) . (34) 

It is thus obvious that there is a semi-spherical structure in ds 2 , which implies some 
form of symmetry of the kinetic energy operator, beyond those associated with rotational 
symmetries. We will soon learn the importance of this symmetry. Moreover, it is easy 
to get the inversion of the transform, namely 

/i = ^(l-£cos/3), /2 = ^(i + ecos/5 ), / 3 = ^ sin/? , (35) 

where £ = sin^. For convenience, we shall use (p, £,/3) instead of (p,9,f3) throughout 
this work, and call them the Jacobi spherical coordinates for reason that will become 
obvious later (see (fTSj) and (JZHJ). 



CONTENTS 



13 



3.2.3. Interparticle distance expressed in Jacobi spherical coordinates From the 
cyclic permutations of three sets of Jacobi coordinate vectors, i.e., 



,(2) 



r (3) 



mi 



1 — mi 



(m 2 (ri -r 2 ) +m 3 (ri -r 3 )), y 



(i) 



1 



m 2 

fflo 



(m 3 (r 2 - r 3 ) + mi(r 2 - ri)) , y 



(2) 



m 3 



y 1 - m 3 
and two simple relations 

( ? ) ■ 

/ x( 3 ) \ 

v(3) " 



mi(r 3 - ri) + m 2 (r 3 - r 2 )) , y 



(3) 



1 m 2 m 3 
1 — mi 

1 m 3 mi 
1 - m 2 

1 m 1 m 2 

1-7713 



^2 - r 3 J 



>3-ri), 



ri - r 2 



cos j3 2 sin /? 2 

— sin (3 2 cos /3 2 

cos /3 3 sin j3 3 

— sin /5 3 cos 



X (D 

yd) 
X (D 



where 



cos /3 2 



cos /3 3 



mim 2 



m 3 + mim 2 
I mim 3 
m 2 + miwi3 



sin/3 2 = - 
sin/3 3 = J 



m 3 



m 3 + mim 2 



m 2 



m 2 + m,im 3 



we can write |rjj| 2 in the Jacobi-spherical coordinates (p,£,fl): 



with 



| r 23 
l r 3l| 
l r 12| 

h = 



hp 2 (I cos p), 
k 2 p 2 (l+£cos(f3 + 2f3 2 )), 
A; 3 p 2 (l+£cos(/3 + 2/3 3 )), 



mi 



m 2 



m 3 



(36) 



(37) 



(3* 



(39) 
(40) 

(41) 



2m 2 m 3 2m 3 mi 2mim 2 

For the special case of helium-like ions with infinite nuclear mass, we take 
_ r x + r 2 _ r 2 - ri 
X "^2~' y "^2~' 
as the Jacobi coordinate vectors, where ri and r 2 are the position vectors of two 
electrons, and the distances become 

1 



l r i 

|r 2 
|ri2 



2 = -p 2 (l-£sin/3), 
2 = lp 2 (l+£sin/3), 



= p 2 (l + ecos/5). (42) 

From (f3HJ) and (|4*2*|) . it is straightforward to express the Coulomb potential in terms of 
(p, £,/3). It will be seen that the Coulomb potential becomes separable in the Jacobi- 
spherical coordinates. 
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3. 3. Symmetries and angular momentum eigenfunctions 

Schrodinger equation is invariant under the spatial rotation and coordinate inversion. 
That means both the total angular momentum operator L 2 (= L 2 + L 2 , + Lf), its 
z-direction component L 3 , and the parity operator commute with the Hamiltonian, 
therefore one has the angular momentum quantum numbers (l,m), where m 
characterizes the azimuthal component of the angular momentum, plus the parity 
quantum number A = 0, 1 for even and odd parities, respectively. In addition, the 
wavefunctions must also satisfy the Fermi-Dirac statistics, manifest as antisymmetry 
under the exchange of two electrons. It should be noted that since the antisymmetry 
applies to the product of spin state with the wavefunction, hence the wavefunction can 
be either symmetric (for antisymmetric spin state) or antisymmetric (for symmetric spin 
state). 

In this section, we use a set of bi-harmonic functions as the eigenfunctions of L 2 
and L3, and expand the wave functions in these basis functions so as to separate out 
the rotational degrees of freedom from the Schrodinger equation. 

The system is independent of the choice of z-direction, hence it is enough to consider 
the special case, m = I. As is well known, y l m (v) = r l Y^{9 , 0) , where the set (r,6,<j)) 
denotes the spherical coordinates of a vector, is a homogeneous polynomial of degree 
I with respect to the components of r, and satisfies the Laplace equation. It is the 
eigenfunction for the angular momentum operator. The polynomials y^ l ('x)y l r ^, q (y) can 
be combined to form the eigenfunctions of the total angular momentum L 2 (L = L x +L y ) 
by the Clebsch-Gordan coefficients jSHj- For the special case, m — I, one has two types 
of bi-homogeneous functions [HI E] with respect to x and y: 



This specific family of basis functions, h a l(x, y) (A = or 1), has the following 
properties. 

(i) They are in the null space of the kinetic energy operator. In particular, they have 
zero eigenvalues for the operators A x , A y and A xy , i.e., they exhibit bi-harmonicity: 




(43) 



where 



u = £1 + ix 2 , v = y 1 + i?/ 2 , 




(44) 




(45) 




dxjdyj 



(ii) They are bi- homogeneity in Xj and yj. 
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(iii) They are the common eigenfunctions of L 2 and L 3 , i.e., 

L 2 hi x i = (a + b-X)(a + b-X + l)hi x l 

L 3 h { a X l = (a + b-\)hi x l (46) 

(iv) They exhibit spatial inversion symmetry under (x, y) — > (— x, — y): 

^(-x,-y) = (-ir + ^(x, y ). (47) 

(v) They exhibit symmetry under the permutation of particle 2 and particle 3, y — > — y: 

^(x,-y) = (-l) b C(x,y). (48) 

By using the ladder operator L_(= Li — iL2), it is easy to construct the common 
eigenfunctons of L 2 , L 3 , and the parity operator from h^: 

C;* = L -C ( 49 ) 
It is simple to verify that the set of functions 

{ht ] q+ X, q ;l-m, A < q < 1} (50) 

are the common eigenfunctions for L 2 , L3, and the parity operator with the eigenvalues 
/(/ + 1), m, and (— l) l+x , respectively. 

3.3.1. Wave function expansion For the special case of m = I, one can expand the 
wave function with parity (— l) l+x in the following form: 

* w = I>?^A>,y), ( 51 ) 

where A = or 1, and [ipg X \ A < q < l\ are (/ + 1 — A)-tuples of functions of the 
rotationally invariant variables (fi, f2, fs). They are the coefficients of the vector space 
|^1j + a,5( x i y)i A < g < /|. This expansion is unique because of the orthogonality of 

j^ +A >, y ), A<g</}- 

3.3.2. Reduced Schrddinger equation These operators of A x , A y , V x and V y , on 
4>{fi, f2, fs) niay be expressed in terms of fi derivatives as 

( d 2 d 2 d 2 d \ 

( d 2 d 2 d 2 d \ 

Aip = (A x + A y ) ip 

f f d 2 d 2 \ ( d 2 d 2 \ d 2 ( d d 

= 4 N + h WV + 4f3 [dKdfs + dm) + (/l + h) df! + 6 [Wi + dh 



v ^= ( 2 ^ + " 1 4' 2y2 4 +s2 4' 2y3 4 +x3 4)^- (52) 
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(A) 



-g+A,g 



I 

E 



A + 4(1 - q + A)^r + Aq^- ) ^ + 2(q - X)^ 1 + 2(Z - q)^ 1 







dfi 



(A) 



dip, 



(A) 



a l-q+\,qi 



(53) 



where A = or 1. Substituting (J53j) into the Schrodinger equation, one obtains a system 
of coupled partial differential equations (PDEs) for the (— 1)' +A parity states, in terms 
of the expansion coefficients: 




A^ + 4(l-q + X)^f- + Aq 



9/i 



9/2 



(54) 



-2(g-A)^-i + 2(/-g):, 

where A = or 1, and A < q < I. For a given angular momentum I , one has Z + l coupled 
PDEs for the (— 1)' parity state, and I coupled PDEs for the (— 1)' +1 parity state. In 
terms of atomic spectrum, these (— 1)' +1 parity states are called doubly excited states 
(DES), such as the spectra 2p 2 3 P e and 2p3d 1,3 D°. Below we note three special cases 
of interest. 

(i) The case of I = 

In this case, the wave function is a rotationally invariant function, and the 
Schrodinger equation becomes 
1 



2M 



Aip + Vip = Eip. 



(ii) The even parity case of I = 1 
In this case, the wave function 



and the coupled PDEs reduce to one single PDE: 



1 f A , M 
{ Aib + 4— - 

2M ) v 



dfx df. 



^) + v.p 



Eip. 



(55) 



(56) 



(57) 



(hi) The odd parity case of I = 1 
In this case, the wave function 

and the the coupled PDEs become 

,dipo 



(5? 



1 

2M 
1 



A^ + 4 



- — < Aipx + 4 
2M I ^ 



Mi 

dijii dipp 
dh df 3 



+ Vip 
+ Vip 1 



Eipo, 
EiPl 



(59) 
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In the above form of (|54j) . the rotational degrees of freedom have been completely 
projected out from the Schrodinger equation. This is achieved by expanding the wave 
function in the bi-harmonic basis, |/i|_ g+A g (x, y), A < q < i|, leading to a set of coupled 

PDEs for the expansion coefficients functions |?/^ A \ A < q < ij, expressed with respect 
to the three rotationally invariant variables (/i,/2, /s)- The problem is hence three 
dimensional. The number of the coupled PDE system is noted to be finite, / + 1 or 
/, and there are no singularities, which is noted to contrast with the case if the Euler 
angles were used. 



3.4- Expansion of reduced Schrodinger equation in J acobi- spherical coordinates 

We introduce the Jacobi-spherical coordinates (p, £, (3) to rewrite the reduced 
Schrodinger equation. It will be seen that the present step naturally leads to the 
introduction of the Jacobi polynomials as the eigenfunctions of the angular differential 
operator, hence the denotation of Jacobi-spherical coordinates. It should be noted 
that the Jacobi-spherical coordinate has been used by Simonov jSH] and Whitten [60J, 
the latter used it to study the expression of pair potential. Mandelzweig and co- 
workers jnUESEH] have also used this coordinate system to calculate some properties of 
the Coulomb system. In the present case, the Jacobi-spherical coordinates are applied 
to the reduced Schrodinger equation for the three-body system. As such, it offers a 
natural coordinate system to delineate the intrinsic symmetry of the reduced three- 
body kinetic energy operator. We first note some physical interpretation and properties 
of this coordinate system. 

(i) The variable p = V% where / = J2^=i mjrj is the moment of inertia in the center- 
of-mass frame, thus p provides a natural measurement of the "size" of the system. 

(ii) For the three sets of Jacobi coordinate vectors of (|36|). there are three simple 
relations among them: 

P, 

(3^ _ 2 /5 3 = p, (60) 

where p and £ are unchanged, and (3 is shifted. This property has been used to 
express r^. 

(iii) For the permutation of particle 2 and particle 3, y — > — y, (p, — > (p, £,—/?). 
This property will be useful for us to deal with the case in which there are identical 
particles in the system. 

In terms of the coordinate system (p, £,/?), one has the following differential relations 
for the function if)(fx, f 2 , / 3 ): 

dip 1 dip 1 
dfi 2p dp p 2 



p (1) 


= P (2) = 


p (3) 




= e {2) = 


£(3) 


/3« 


= (3^ _ 


2/3 2 



dip 



dip sin (3 dip 



cos/3--— — 
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and 



where 



dip 

W2 
dip 



1 dip 
2p"dp 



dip ( dip 



— — = — sin/3— + 



P 1 



dip cos /3 dip " 



sin /3 $0 " 



p- 



d£ 



a/r p dp p z 



3£ 2 ^ 1 d 2 ip 



(61) 



(62) 



(63) 



$.^..Z. Reduced Schrddinger equations in terms of (p, £,/3) In general, there are two 
quantum states of different parities for the special case of m = I, whose wave 
functions are given by = Y,q=x^n^h\_ q+Xq (y^,"y), where A = or 1. Here 

Up^\ A < q < /} are (/ + 1 — A)-tuples of rotationally invariant functions satisfying 
the systems of coupled PDEs. Substitution of the derivative relations, (JoTj) and (JfSJ), 
into the reduced Schrddinger equation, the coupled PDEs of (J53|) . leads to the following 
coupled PDEs in the Jacobi- spherical coordinates (p, £,/?): 



1 

2M 



r d 2 



dp 2 



+ 



5 + 2/ + 2A dip^ 



P 



dp 



+ 



P 1 



d 2 ip^ 1 
i — — | — 

de 



(l + X + 3)edip { q x) 



+ -^r- + (l-2q+ X)B 1 iP^ + {q- A)£ 2 ^\ + (I - q)B 2 ip { q x } l 



ViP, 



(A) 



£ 2 d/3 2 

(A = 0, 1;A < q < I) (64) 

The above system of PDEs has remarkable simplicity and uniformly. Notice that the 
differential operators involving partial derivatives in p are the same for each component 



d 2 4 X) , 5+21+2X d4 X) . 
'q ) > dp 2 "T" p 9p ' 

organized into the following two parts. 



function i.e.. 



while the partial derivatives in £ and /3 can be 



(i) A uniform part for each component function ip q x \ denoted A^ l+X \ which comes from 
the angular part of the kinetic energy operator (see 



d 2 



1 - (/ + A + 3)£ 2 d 1 <9 2 



[ii) A second part consisting of (I — 2q + X)Biip^ + (q — \)B 2 ip q X ? 1 4- 
where 

9 cos/3 d 



d sin/3 9 

-Di = COSD— — — , 

d£ £ dp' 



B 2 = sin (3- 



(65) 

(I - q)B 2 iP q % 
(66) 



di i dp 

This part comes from the first order derivative terms (with respect to f\ ) in (|54|). 

Below we write out the reduced Schrddinger equations explicitly for those cases of 
particular interest. 



CONTENTS 
(i) The case of / = 0: 



1 (d 2 9 5<9* 4 



2M [dp 2 pdp p 2 
(ii) The case of I = 1 and of even parity: 



±1*9 9?9 4 (2) ' 

2M 1 dp 2 p 9p p 2 



y* = 



V* = E9. 
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(67) 



(6* 



(in) The case of / = 1 and of odd parity: 



1 \d 2 ^ 7djj 4 
2M \ dp 2 p dp p 2 

_l_(d± 7c^i _£ 
2M {"dp^ + ~p~~dp + p 2 



A^o + + B^A \ + V^ = E^ 



- £1^1 + 5 2 ^ol ^ + Vipi = Eijjx. (69) 



3.4-2. Eigenf unctions of the angular differential operator A® Here we give the 
eigenfunctions for the following type of angular differential operator: 



A iH = (1 _ e 2 )^ + ±±l3^°± + 1^ 



(70) 



Let P^ a ' b \x) be the Jacobi polynomials, then the following set of functions constitutes 
a complete family of eigenfunctions of A®: 

'j { ±n{^P) = e ±i ^r^' m) (2e 2 - 1), m,n > o) , (71) 



(72) 



with their respective eigenvalues given by 



4n ^1 + ^ + m + nj + m(l + m + 2) | . 

(0 



The details and some properties of J± mtn are given in Appendix A. 

3-4-3. Matrices B\ and E>2 applied to the eigenfunctions of A® Applying B\ and B2 
on the eigenfunctions of A^ l \ {j±m.,n(£> (3), m,n> o|, one obtains the following three 
cases. 

(i) The case m = 0: 

n 

r r(0 _ 9 >T fcW / rW 1 r(0 \ 



fe=l 



i R T(0 _ O V /.W J f W f W 1 

1 - D 2»'o,n — z °0,n.,fc \ J l,n-fc J -l,n-k] • 
k=l 

(ii) The case of positive indices: 

/ n n 

r r(0 _ 9 [ V /> W 4- V h W T {1) 

±J ^- U m,n * \ Z_^i u m,n,k u m+l,n~k ' u m,n,k u m-l,n- 



(73) 



k 1 



\k=l 
' n 



j r t(0 =2 fV b {l) T (l) - V b [l} T [l) 

i J->2< J m,n I / , m,n,k m+l,w— fc / t m,n,k m—l,n 



k=0 

n 



(0 r(0 



(74) 



fc=0 
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(iii) The case of negative indices: 

/ n n 

-"IV-m.n — * \ 2^, u m,n,k J -(m+l),n-k ~T U m,n,k J -(m-l),n-fc 
fc=l fc=0 > 



(n n \ 

V-l(0 t(0 /- (/) 1 f7^ 

°m,n,fc J -(m+l),n-fc °m,n,k J -(m-l),n~k j ■ K<°) 

k=l k=0 I 

In the above, the coefficients fc and b® n k are given by 

bm, n ,i = ^{a + m)+n, 

.(;) m + 71 / 1 . . 
brn,n,2 = ~ , -(a + m) + 71 - 1 

a + m + n \2 
a) m + n m + n — 1 (1 

b m,n,3 = — : 7 ~ 7 \ 7<( a + m ) + n ~ 2 > ete - 

o + m + no + m + n-1 \2 



~m m + n ( 1 



m + n m + n — 1 /l 

a + m + na + m + fi-1 V2 
m + n m + n — 1 m + n — 2 /l 



& m,n,i = — : : : : ~ \-{a + m)+n-l 

SO 



C.9 = — ; ; ; : 7 — ; ; x \ -(a + m) +n-2 \ , etc. 



m ' n ' 2 a -|_ m _|_ na -|- m -|-n — la + m + n — 2 \2 

(76) 

where a = ~. The details of calculation can be found in Appendix A. Denote 

|m ) n;/> = Cg B jg B) (mGZ,nGZ+), (77) 

where is the normalization constant of J$ n , Sf 1 and being the matrices of 
operators B\ and _B 2 on the function space {\m, n\l) ,m G Z,n G *6i and $2 

have two essential properties: 

(i) For the index m, both &i and B<p are tridiagonal block matrices with zero diagonal 
blocks. 

(ii) For a fixed m, each none-zero sub-block matrix of S x and for the index n is 
an upper triangular matrix. 

The above two properties will be used to advantage in designing the numerical scheme. 

3.4-4- Coulomb potential expressed in J acobi- spherical coordinates (p, £, (3) From (|19|) 
and (jnnj), it is easy to get the potential in the Jacobi-spherical coordinates 

V=^Ml, (78) 

where 

1/2 /„ X 1/2 /„ x 1/2 



2 2 { 277127723 ^ Z Z { 2m 3 m i ^ Z Z { 2mimg 

rj(e, /?) = 2 H 1 "™ 1 4 + ; 3 1Vl " m2; + ; 1 2 ^-"^ (79) 
Vl + ^cos/^ /1 + £ cqs^ + 2 /3 2 ) ,/l + £ cos(/3 + 2/5 3 ) 



From (fT8|) and (|79|. one can see that 
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(i) The Coulomb potential is separable into p and (£, (3) components, which is crucial 
for further reduction. 

(ii) The potential has three singularities in the (£,/3) space: (1, 7r), (l,7r — 2/3 2 ), and 
(1, 7T — 2j3 3 ). These singularities affect the convergence of the numerical calculations. 

3.4-5. Matrix elements ofU(£,/3) with respect to the eigenfunctions of A® The matrix 
element of [/(£, (3) in the space {|m, n; I) , m G Z,n G Z + } is defined as 

ZY« (mn, m'ri) = f C exp (-i(m - m')/W6 /3)T£„ m , n , (£) d£d/3, (80) 
where 

r£, mV (0 = Cg B C7SU W+|m ' l+1 (l -0 ,/a ^ lAH) (2^ a - l)Pr' |m ' D (2f - 1). (81) 
Setting 

?7 (e,/3) = (l + ecos/3)- 1 / 2 , (82) 

and denoting its matrix Uq (mn,m'n') in the space {\m,n; I) ,m £ Z,n E Z + } the 
"universal matrix" , one has 

U^\mn,m'n') = C(m — m')UQ l \mn,m'n'), (83) 

where 



C(n) = J^HUZA + tf^Z& exp(,2„ A ) + . ^Z,Z 2 exp(i2 nA ) (84) 
V 1 — mi VI - m 2 VI - m 3 

is denoted the "modulation factor" . A few properties of the matrix is noted below. 

(i) The physical quantities, i.e., masses and charges, appear only in the modulation 
factor. For the special case of particle 2 and particle 3 being identical (rri2 = = 
m , Z 2 = Z 3 = Z and (3 2 = ~ P3 = Po), the modulation factor is a real number, 
given by 



C(n) = Jm~ Z 2 + 2J^^ZZ 1 cos(2n/3 ). (85) 
V 1 - m 

(ii) The universal matrix is independent of any specific system, which is the nature of 
three-body Coulomb system. In the numerical scheme, this property will be used 
to advantage. 

In the particular case of helium-like ions with infinite nuclear mass, we would like to 
note that the Hamiltonian is 

H = -i(Ai + A 2 ) + (-!L_Z + ±_\ (86) 
2 V r 1 r 2 r 12 J 

where Z is the nuclear charge. In this case the modulation factor becomes 

C{n) = 1 - 2v / 2^cos(nvr/2). (87) 
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3. 5. Further reduction to a system of ID ODEs 

The "space-rotation reduction" enables us to completely remove the three rotational 
degrees of freedom, and reduce the Schrodinger equation to a system of PDEs solely 
in terms of the rotationally invariant parameters (/i,/2,/s)- In terms of the Jacobi- 
spherical coordinates (p, £,/?), the above reduced system of PDEs becomes well- 
organized, in the sense that there can again be a separation of the angular component 
with the radial component. The angular component consists of the sum of A^ l+X ^ and 
the integral multiples B 1 and B 2 . The complete solutions of the eigen-f unctions of A® 
in terms of the Jacobi polynomials, namely |j±m,n(C> /?)}> can t» e explicitly obtained, 
as well as analytic formulas of the matrices of the linear differential operators B\ and 
B 2 with respect to the above eigen-basis of A"'. Furthermore, we provide a way of 
computing the matrix of U(£,/3), which constitutes a separable part of the Coulomb 
potential. By combining all the above results, it is straightforward to further reduce the 
system of PDEs in (p, £, j3) to an infinite system of ODEs in terms of p. By writing 



E E4S»K W ^ + A ): (A = 0,l;A<g<0 



-oo n=0 



and substituting the above into the system of PDEs (|64)l . we obtain a set of ODEs on 
the expansion coefficients. That is, let !F^ x \p) be the column vector with {/^.^(p)} as 
the components, and let B[ l+X \ B 2 , U^ l+X ^ be respectively the matrices of the linear 
operators B\, B 2 and U(£,j3) with respect to the eigenfunctions of A^ l+X \ Then the set 
{j~q X Hp)i A < g < /} satisfies the following system of coupled ODEs: 



1 

2M 



d 2 , 5+21+2A d \ WA) 
dp 2 J p dp) 1 
~ ^(*+A)jr(A) 



{I- 



2q + \)B[ 1+X) FW 



\)B { 2 l+x) T {x \ 



(i+A) x-(A) 
9+1 



p 



(89) 



where A = or 1, A < g < Z, and *4.^ +A ) is the diagonal matrix of the eigenvalues of 
A" +X K In particular, we note the following special cases. 

(i) The case of / = 0: 




fii) The case of I 



fifi) The case of I 



2M 



2M 



a 2 

dp 2 

+ p 2 

dp 2 
+ ^ 



1 and even parity: 
1 U d 2 9d_ 
2M \ [dp 1 + pdp 
1 and odd parity: 




-U^T = ET. 
P 



t + - 2 a^t\ 



-U^T 

P 



ET. 



(90) 



(91) 



pdp 



A^T + B\ 1} T 



^(i)^ _ b^T 



(92) 
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By using the following notations 







K ) 




B? 









B= ' A _,„ , U 





o 


O 


AM 




O 


O 





(93) 




ET. (94) 




where O represents the zero matrix, and A, B and U are all block matrices, the 
two coupled equations (fQ2j) can be abbreviated as 

--((- 

2M \ \dp 2 

In the general case of / > 2, we can use the block-matrix notations to write the system 
of ODEs in a more compact form: 

where A = or 1, JF^I is the collection is the collection of 

the matrices A^ l+X \ B^ is the collection of the matrices Sf +A ' ) and S2 , and is 
the collection of the matrices U {l+X) . All of A [x] , B [x] and U [x] are (/ + 1 - A) x (/ + 1 - A) 
block matrices. As it turns out, these block- matrix notations are not only convenient 
for analyzing the dependence of J-'^(p) on p, but they also provide more insight into the 
structure of the three-body system, which will greatly benefit the numerical calculations. 

At this point it should be noted that while the starting point of the three-body 
problem is a nine-dimensional problem, in ()95j) the problem has been reduced to a one- 
dimensional problem, since J 7 ^ (p) depends only on p. Below we provide a way to solve 
this one-dimensional coupled ODEs by converting them into a linear eigenvalue problem. 

3.6. Conversion to a linear algebraic eigenvalue problem 

In ()95|) . there is a common differential operator, + 5+2 ^ +2A A ; f or each component of 
the vector T^ x \p). The energy level E of a stationary state should be negative, and if 
we analyze the "asymptotic behavior" of (|95j) . the limiting equations are all reduced to 
the following single ODE: 

hw- B )^° (£<0) - (96) 

Equation (J96)) clearly implies the exponential asymptotic decay of the wave functions of 
stationary states and, moreover, the order of such exponential decay is given by 



k = V-2ME, E < 0. (97) 
In other words, 

F [x] ~ exp(-fcp), as p -> oo. (98) 
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Setting x = 2kp, F [x] (p) = e~ kp y^(x) and 

' d 2 21 + 5 d 



then 



pe kp L {l+x \p) = 2k 



dp 2 p dp 
d 2 



x 



+ (a + 2A + 5-,)jL-(/ + A + §) bw 



(99) 



(100) 



3.6.1. Expansion of the ODE system in terms of the Laguerre polynomials. If we 
substitute the above results into (j95j) . the system is transformed into a system of ODEs 
with respect to the new independent variable x: 



where A = or 1, and 
Id) 



y 



A] 



4 



x 



A [x] y [x] + B [x] y [x] 



M 
~k 



d 2 x d 

+ (2/ + 5 - a:) — . 



(101) 



(102) 



da; 2 ' v ' dx 

These Laguerre polynomials ^L p 2l+4 \x), p = 0, 1, 2, • • -j are the eigenf unctions of 



i.e., 



,r- 



+ (21 + 5 - x) 



d ] L (2«+4) 



-p 



r (2*+4) 



(103) 



y dx 2 ' v dx J p 

Equation (|103j) shows that it is advantageous to use the family of Laguerre polynomials 
[Lf l+2X+A \x), p = 0, 1, 2, • • •} as the basis functions to solve (fTTTTJ) . 

3.6.2. Infinite linear equations for the coefficients. Let be the undetermined 
coefficients in the expression of y^(x) as the linear combination of the Laguerre 
polynomials |l^ 2Z+2A+4 ^(x), p — 0, 1, 2, • • - j: 

oo 

j;W( a ;)=^ u WLf +2A+4 )(x). (104) 

p=0 

Here w£ A ) is more precisely defined as the set ^v^ m n , A < q < iX. By omitting the 
sub-indices q, m and n in the block-matrix notations, we obtain 

r -E P °lo(p + ' + A + f)4 A) 4 2i+2A+4) ^ 



L A [X] 



- M V ?> (A) /-( 2i + 2A + 4 )7yW 



(105) 



where A = or 1, and Z^a] is the identity matrix, with the same rank as A^. 
Multiplying the above equations by kx and then making use of the recurrence relations 
of the Laguerre polynomials, we get the systems of linear relations in terms of the 
Laguerre polynomials. By applying the orthogonality relations of these polynomials, a 
system of algebraic linear equations for the coefficients {i^^^j is obtained. In the 
block-matrix notations, they may be expressed as 



k 



p(p + i + A+|)^ A _ ) 1 
-2 (p + I + A + | 



(A) 



+ (p + I + A + |) (p + 21 + 2A + 5) 



/A) 
} P+1 









> = m| 







(A) 



+ (2p + 21 + 2A + 5) v p x) }U [X] (106) 



■ (p + 21 + 2A + 5) v 



(A) 
P+l 
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1) 






T> ( p\p,p + 1) 


Gp\p,p- 


l) 


Gp(p,p) 




Gp(p,P + 


l) 



where A = or 1. Let T^p and Qp^ be the tridiagonal matrices related to the index p, 
and Tp) be the identity matrix with the same size of T>p and Qp\ where 

P (p + 1 + 2 

p + l + ^ (p + 2Z + 5) . 
(2p + 2/ + 5), 

-(p + 2/ + 5). (107) 
Equation (|1U6|) can be re-written into a compact form: 

4S ? '«1> + ««)}" < ^ M ^ +A)8M[X '^ <A ' (108) 

where the symbol ® denotes outer-product between two matrices. Three points should 
be noted. 

(i) Equation (|108|) has a robust structure, mainly due to the fact that p and (£,/?) are 
separable in the potential. Vp +X ^ and C/p +A ^ are unchanged when there are identical 
particles in the system. 

(ii) The left-hand of ()108j) can be partitioned into two parts, namely 

M = V [ p +X) ® T AW + 4J? +A) ® A [x] , Mi = 4X { p +x) ® B [x \ (109) 

where 7W0 is a super-tridiagonal matrix that dominates the left. Consequently, we 
can solve the linear equation, M. Q x = y, just like solving a tridiagonal equation. 
Furthermore, the equation (A4 +Aii)x = y can also be solved based on an iteration 
process, detailed in the next section. 

(iii) The outer-product form provides a way to easily and efficiently get the matrix-on- 
vector product. 



4. Numerical scheme 

Schrodinger equation for the Coulomb three-body system has been reduced from a nine 
dimensional problem to a one dimensional problem and finally to a linear algebraic 
problem in section El Here we present a numerical scheme designed to solve the relevant 
linear equations. The core of the scheme is the iterative solution of the matrix-eigenvalue 
problem. The advantage of the iteration procedure is that it only requires the matrix- 
vector product, which can significantly facilitate the numerical calculations due to the 
sparse-block structure of the matrices in our problem. The general linear eigenvalue 
problem is written into a standard form, and by applying a special integration rule, 
presented in Appendix B, the calculation of the potential energy (Coulomb interaction) 
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matrix elements and the product of that matrix with a vector are combined together, to 
achieve a high degree of numerical accuracy with minimal computational resources. In 
present numerical calculations, the aim is to solve a system of sparse linear equations, 
and an iteration solver is adopted for that purpose. Since truncation is necessary, an 
extrapolation procedure is applied to accelerate the convergence. 

4-1- Solution of the linear eigenvalue problem 

After introducing the Jacobi polynomials, the system of PDEs is reduced to a system of 
ODEs, and the ODE system is further reduced to a linear eigenvalue problem through 
expansion in terms of the Laguerre polynomials. In the block-matrix notations, the 
linear algebraic equation can be expressed as 

k {Vl +X) ® Z m + 4J? +A) ® (aW+BW)}vW = M{S? +A) ®Z/M}vW (110) 

where the symbols' definitions can be found in section El and k is related to the energy 
Eby 



k = V-2ME, (111) 

where M is the total mass of three particles. Equation (jllOj) may be compactly expressed 

as 

k M M L v = M R v, (112) 
where ku = k/M, and 

M L = Vf x) ® 1 AW + 4J? +A) ® (AW + £ [A] ) , 

Mr = g { p +X) ®U [X] . (113) 

In general, only the low-lying states are of interest, which implies only a few of the 
largest eigenvalues of ()112|) are needed. One can write this general eigenvalue problem 
in a standard form: 

Ml l M R v = k M v, (114) 

where Ai^ 1 means the inverse of Ml- Hence, the task becomes to design a scheme to 
find a few of the largest eigenvalues of l}114j) . 

We use one of the software packages, such as ARPACK jM], designed to compute a 
few eigenvalues and their corresponding eigenvectors of a general n-by-n matrix A. It is 
most appropriate for large sparse or structured matrices where structured means that a 
matrix- vector product w <— Av requires order n, rather than the usual order n 2 , floating 
point operations. This particular software routine is based on an algorithmic variant 
of the Arnoldi process called the Implicitly Restarted Arnoldi Method. By using this 
package, only the operation of the matrix-vector product is needed. In our problem, the 
following two steps are used to calculate the matrix-vector product M^M-rv. 

(i) First calculate y <— M.rv. This step is the most time-consuming part, and the 
block- matrix structure of Mr is found to be suitable for distributed computing 
systems, such as PC clusters. 
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(ii) Calculate y <— M. L l y. The difficult part of this step is to solve a large sparse 
equation, and a sparse solver is adopted. 

Below we detail each of the two steps. 

4-1.1. Matrix-vector product (M.rv) evaluation The matrix Mr, derived from the 
Coulomb potential, is the tensor formed from the outer product of Qp +X ^ and IA^: 

M R = gp X) ®U [X \ (115) 

where Qp +X ^ is a tridiagonal matrix, and IA^ is a block-diagonal matrix in which every 
sub-block is the same, denoted U^ l+X \ 

The Kronecker tensor product, X cg> Y, of two matrices is a larger matrix formed 
from all possible products of the elements of X with those of Y. If X is m-by-n and Y is 
p-by-q, then X ® Y is an mp-by-nq matrix. The elements are arranged in the following 
order: 

/ * Y X 1}2 * Y ... X 1>n * Y \ 
X 2 a*Y Xo 2 *Y ... X 2 „*Y 



X®Y 



(116) 



\ X mA * Y X m>2 * Y ... X m , n *Y J 

For the matrix-vector product {X ® ^}v, where v is a vector, it is advantageous to 
reshape v to a q-by-n matrix V, and to calculate R = (Y *V) *X T . The vector reshaped 
from the p-by-m matrix R is the final result. 

The particular structure of A4r and IA^ makes it easy to calculate the matrix- 
vector product. The tensor form of M.r provides a straightforward way to calculate 
the matrix-vector product, and the block structure of IA^ lets one to focus on only the 
sub-block matrix IA^ 1+X \ The elements of matrix lA^ are expressible as 

IA^ l \mri,m'n') = C(m — m')IAQ l \mn,m'n'), (117) 

where both the "modulation factor" C{n) and IAq the "universal matrix" are given 
previously. From the definition 

" x cos (m — m')(3 

//3=q/^=0 y/1 + £cOS/3 

where 

Mt) = £(1 - ef\ P®(£) = CS n £ H ^ /2 ' H) (2£ 2 - 1), (119) 
with jp^/ 2 'l m l)(x) j being the set of Jacobi polynomials. Set 



lA«\mn, m'n') = 2 f f 008 ( ™ ^ w^) P£ n (0 ^U) d£d/3, (118) 
J3=oJf=o v 1 + c cos p 



and 



«m(0 = f cos m/3W 1 ^ d/5, (120) 

Jo V 1 + £ cos p 



Dff({) = 2C ' m ^g>"-«) , (1211 
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then (|117|) becomes 

UU(mn,m'n') = f 17^,(0 p£(0 P$n'(0 d£. (122) 

The integration implied by ()122|) is performed using the IMT integration scheme (see 
Appendix B). 

4-1-2. Calculation of Ai^v One needs to solve a linear algebraic equation to get the 
vector AAj^v. From the previous discussion, it is known that AAl is composed of two 
parts: 

M L = M + M 1 , (123) 

where M = V { p +X) ®l AW +4Zp X) is a super-tridiagonal matrix that dominates 

the matrix j\4l, and Ai\ = 4Xp +A ^ <g> is relatively small. Because of the super- 
tridiagonal structure of Aio, one can obtain the vector A4q V exactly like solving a 
tridiagonal equation. Obtaining Ai^v may be implemented as 

Ml l w = (J + M^M^M^v, (124) 

where I represents the identity matrix, through the following two steps. 

(i) First solve for 

y <— Mv V (125) 
As described above, this step can be easily implemented. 

(ii) Solve for 

y^(l + M 1 M 1 )- l y. (126) 

Since the matrix Aio dominates, thus A4q l Ai\ is small in some sense, so one can 
consider the matrix X + M.q 1 J\Ai to be close to the identity matrix, i.e., 

l + M^Mx ~ J + e. (127) 

Thus an iteration solver should be efficient. In our calculations, the Conjugate- 
Gradients-Squared (CGS) jnS] method is taken as the solver. The CGS method 
just requires the user to provide the matrix- vector product, i.e., one works out the 
vector 

(l + Mo l Mx)v (128) 

during the intermediate iteration steps by first calculating t <— A4±v, and then 
calculating t <— A4q H. The final result is given by v + t. In practice, after about 
12 iterations, the relative error of the solution is reduced to 10~ 14 . 

Since the infinite matrices in the problem are truncated into finite matrices, the 
truncation is done by following a rule, specified below. Results obtained from smaller 
scale calculations, i.e., the eigenvalues and their eigenvectors, are taken as the initial 
guesses for the subsequent, increasingly larger cases. This process fully takes the 
advantage of the previous results, and significantly speeds up the calculation. It also 
provides a sequence of data from which one can apply the extrapolation procedure (see 
below) to obtain a few more digits of accuracy. 



CONTENTS 



29 



4-2. Truncation and extrapolation procedures 

In our numerical calculations, the expansion of the partial wave function ipw has to be 
truncated, i.e., 

M N P 

$ A) = e-* E EE^l^S A He,/5)4 2 ' +2A+4) (2M, (129) 

m=—M n=0p=0 

where M, N and P represent the three truncation indices for the m, n and p, respectively. 
One solves the finite dimensional algebraic eigenvalue equation to obtain the lowest order 
eigenvalues and their corresponding eigenvectors. These eigenvalues are functions of M, 
N and P, i.e., 

E — E(M, N, P), (130) 
where the desired answer is the limiting value 

E^ = E(oo, 00, 00). (131) 

4-2.1. Truncation path We have carefully checked the dependence of E on the three 
numbers M, N and P. In our procedure, we chose a "path" in the {M, N, P} space 
to approach the limit E^. First, it was found that for sufficiently large P (typically 
P ~ 70), increase in P does not provide any significant improvement on E. This 
behavior owes to the fact that P controls how many terms in the Laguerre polynomials 
are included in the expansion, but there is a direct relation between the spatial extend 
of the system and the spherical variable p, so for large enough P the terms of Laguerre 
polynomials in the expansion are sufficient to delineate the spatial domain. Different 
states have different "appropriate" P values, and after that value is determined, denoted 
Po, the following "path" in {M, N} space is taken: 

M = M + 20K, N = N + 10K, K = 0, 1, • • ■ , (132) 

where Mq and No are the initial truncations for m and n, generally in the range of 70 
and 50, respectively. The three-dimensional sequence E(M, N, P) is thus mapped to a 
one-dimensional sequence: 

E{K) = E(M + 20K, N + 10K, P ). (133) 

There is a reason why we increase M twice as fast as N. It is known that the eigenvalues 
of operator have the form 

Ag n = - |4n ^1 + l - + m + rij + m(l + m + 2) j . (134) 

As m and n approach infinity, \$ n has the asymptotic form 

A« n oc-(m + 2n) 2 , (135) 
so that the stipulated sequence is designed to coincide with this asymptotic behavior. 
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4-2.2. Extrapolation procedure For extrapolation, an auxiliary variable x, defined as 

x 100 - 100 (136) 

M + 2N M + 2iVo + 40ir v ; 

is used to determine the form of E(x) as provided by the data sequence. A power law 
form is assumed, such that 

E(x) = E 00 - c x a (l + Cl x + c 2 x 2 + ••■)• (137) 

The task is to identify the leading term x a . Sequences of E(k) (E(x)) are collected in the 
numerical calculations. The leading term is identified by studying the relation of E(x) 
versus x at , where at is an estimated value of a. If the plot is a straight line, the leading 
term can be immediately identified. The process of Richardson extrapolation [HE] is 
then performed on the sequence {E(xi),i = 0, 1, • • •} to estimate the limiting value E^. 
The Richardson extrapolation utilizes the following sequences: 

E{x ) = E OQ - c x% - Cl x a +1 - c 2 x a + 2 

E( Xl ) = £oo - coxl - cix? +1 - c 2 x a + 2 

E(x 2 ) = E OQ - c x a 2 - Cl x a 2 +1 - c 2 x a 2 +2 

E{x n ) = E oa - c < - c lX a n +1 - c 2 x a n +2 , (138) 

to approach E^. By defining 

gi {n) = x a n , g 2 (n) = x a n + \ g 3 (n) = <+ 2 , • • • , (139) 
and denoting E( E„, (I138D becomes 

E n = E 00 - c gi{n) - Cxg 2 {n) - c 2 g 3 (n) , n = 0, 1, 2, • • • . (140) 

Then the recursive E — algorithm declares that the limiting value -E^ can be approached 
by the following rule. Let 

E^ = E n , n = 0,l,--- 

9$ = 9i(n), n = 0, 1, • • • and* = 1, 2, • • • . (141) 
For k = 1, 2, • • • and n = 0, 1, • • • one has 

fit" = - 1£) _ w ' ■«&■>■ < 142 > 

9k-l,k 9k-l,k 

where the g^Z 1 fc 's are auxiliary quantities recursively computed by 

9 ( k1 = 9 { k% ~ "til ~ 9 Z'' k ■ 9 ( k%k, * = k + l,k + 2, ■ ■ • . (143) 

9k-l,k 9k-l,k 

The sequences {E^\n = 0, 1, • • • , and k > 1} are more convergent than the initial 
sequence {Eq ,n = 0, 1, • • •}. From the convergence pattern of {E^}, one may 
determine the limiting value with high accuracy. 




Figure 1. The triangle formed by three particles. 



5. Results on three-body systems 

In this section we present results on some three-body systems, and discuss the properties 
of their relevant wave functions. In order to display the data in physical space, we first 
give a description of the relevant angles, as well as the procedure by which the wave 
functions are re-constituted from numerical data. 

5.1. Euler angles and wave functions 

5.1.1. Euler angles and the distribution function The three-body wave function 
depends upon the shape of the triangle formed by three particles and the three Euler 
angles jHHl, such that 

tf = #(A;fi), (144) 

where A represents the triangle which can be described by the three spherical variables 
(p, £, P), or its three edges (r±, r 2 , r 3 ), or (r 2 , r 3 , cos#) as shown in figure[TJ and fl denotes 
the Euler angles that depict the rotational orientation of the triangle in space. Here 
r2 and r 3 denote the distances of particles 2 and 3 from particle 1, respectively, and r\ 
denotes the distance between the two identical particles 2 and 3. The cos# is between 
particles 2 and 3 (the identical pair). For a six-dimensional function, visualization is 
an issue. The Euler angles are noted to take part in the wave function through the 
angular momentum eigenfunctions, but do not enter in the rotation-invariant functions, 
which indicates that in order to deal with the Euler angles, we just need to focus on the 
angular momentum eigenfunctions ? (x, y), A < q < l\. After the Euler angles 

are integrated out from the wave function, a three-variable probability distribution 
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function is obtained, i.e., 

P(A) = / |$(A;fi)| 2 dfi, (145) 
Jn 



which can be easier to visualize. Let P(A) = g{r\, r 2 , r 3 ). In order to display 
the information contained in g(r\, r 2 , r 3 ), we propose to use the following probability 
distribution functions (PDF). 

(i) One- variable PDF, e.g., the distribution of r\: 

P(ri)dri = / / g(r 1 ,r 2 ,r 3 )dr, (146) 

Jr 2 Jrz 

where dr is the volume element. 

(ii) Two- variable PDF, e.g., the radial correlation function: 

P(r 2 ,r 3 )dr 2 dr 3 = / g(n,r 2 ,r 3 )dr, (147) 

Jri 

or the conditional distribution 

P(r 2 ,r 3 ;ri = a)dr 2 dr 3 = g(n = a,r 2 ,r 3 )dr/dr 1 , (148) 

where a is a given value, and a PDF is obtained for each given a. 

We noted that the volume element dxdy in different coordinate systems has the following 
explicit form: 

dxdy oc p 5 £dpd£d/3df2 oc r 2 rldr 2 dr 3 dcos9dQ = r 1 r 2 r 3 dr 1 dr 2 dr 3 dr2. (149) 

In calculations, we always normalize the wave function, i.e., 

f\V\ 2 p 5 £dpd£d/3dn = 1. (150) 

If there are two identical particles in the three-body system, we identify them as particles 
2 and 3. 

5.2. Negative hydrogen ion 

The negative hydrogen ion H~ is an interesting special case of helium-like ions (Z=l). 
It is marginally stable against dissociation into a neutral hydrogen atom plus a free 
electron. The dissociation energy J of H~ ground state is only about Q.lbeV and this 
ion possesses no other bound state [H7| 

The negative hydrogen ion has been found to be of great importance for the opacity 
of sun's atmosphere. The ionization potential of H~, J m 0.75eV, corresponding to 
about 8700° K, is only slightly higher than the temperature for the solar atmosphere. 
As free electrons are released by the ionization of the metal elements present in the 
gas, and since neutral hydrogen is by far the main constituent of the solar atmosphere, 
many of these electrons will be captured to form an abundant source of H~. The 
radiation flux coming from the sun's interior would be absorbed by the H~ ions, 
accompanied by their dissociation. The electrons released can again be captured by 
H~ atoms with the emission of radiation, and so on. The process H~ # H + e~ is the 




r 



Figure 2. Distributions of t\,vi ( or 7-3) and cos(? of the ground state of °°H~. The 
solid line is the distribution of r-i or r% (same), and the dashed line is the distribution 
of r\. Inset is the distribution of cos0. 

main source of observed opacity in the solar atmosphere. The absorption coefficient of 
H has been studied extensively and used in the theory of the solar atmosphere. In 
fact, discrepancies between early calculations and observational evidence on the sun's 
radiation have pointed out the inaccuracies of the calculated H~ wave functions available 
then. Calculation of the wave function and energy of H~ is also of purely methodological 
interest, since this most loosely-bound of all Helium-like ions provides a severe test for 
the various approximation schemes. 

5.2.1. The ground state Bethe was the first to give an unambiguous proof of this ion as 
a bound system (HHl- Using the Hylleraas variational wave functions, Bethe concluded 
that the resulting Rayleigh-Ritz upper bound on the energy lies below -0.5 a.u.. In 1944, 
Chandrasekhar introduced a two-parameter trial wave function, 

\[> = exp(-ar 2 - (3r 3 ) + exp(-ar 2 - (3r 3 ), (151) 

and showed that the energy minimum at a = 1.03925 and (3 = 0.28309 is sufficient 
to provide binding for H~. The function shown in (j!51j) exhibits the specific nature of 
electron-electron correlation in the ground state. The two electrons are on very different 
footings, one bound much closer to the nucleus than the other, which is weakly held 
at a distance ~ 4 — 5 from the nucleus, and this electron can be regarded as weakly 
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Figure 3. Radial density probability distribution P(r) for °°H ( 3 S), H(ls) and H(2s). 

bound in a short-range attractive potential well. In modern variational calculations, 
many r\ terms are necessary in the trial wave function in order to fully account for the 
electron-electron correlation. While many-parameter variational calculations can give 
great accuracy for the ground state energy, the best experimental values come from a 
high resolution (0.03cm" 1 ) laboratory photodetachment laser experiment. The binding 
energy has been determined to be 6082.99 ± 0.15cm _1 for H~ and 6086.2 ± 0.6cm _1 for 
the similar D states [70]. In tabled our results on the ground states of °°H _ and H~ are 
compared with other calculations. In general, all our results have 10 ~ 12 significant 
figures when double precision programming is used. The accuracy in this case is less 
than those obtained variationally. It is well known that the precision of the operator 
expectation values in the variational calculations usually has two less significant figures 
than that for the energy. Thus much more efforts is required to obtain the wave function 
expectation values to the same accuracy. In our case the wave function is calculated at 
the same time as the eigenvalues. No extra effort is required. 

Figure |2] shows the distributions of r\ , r 2 , and cos 9 of the ground state of 
°°H _ . The distribution of cos 9 is a curve that shows a maximum at 9 = n. If the shell 
model is valid, the curve of the cos 9 distribution should be close to a horizontal line 
at 0.5. For comparison, we also plot the radial density distribution P(r), for °°H _ ( 3 S), 
H(ls) and H(2s) in figure El A long tail is seen in the distribution of r 1; r 2 (or r 3 ). In 
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Figure 4. Distributions of n, r 2 (or r3) and cos 6* of the 3 P e state of °°H _ . The 
solid line is the distribution of r 2 or r 3 (same), and the dashed line is the distribution 
of 7*1. Inset is the distribution of cos (9. 

figure [3J the electron distribution of °°H _ is seen to be very different from the neutral 
Is or 2s state. Thus the electron-electron correlation plays an important role in H~, 
and independent-electron model can not give an accurate picture. 

5.2.2. The 3 P e state The even-parity 3 P e state of °°H~ and H~ is quasi-stable, and 
its energy is just below the n=2 threshold (-0.125 a.u.) of the hydrogen atom. For 
the °°H _ ion, the existence of the even-parity 3 P e state was predicted computationally 
nearly 40 years ago, followed by many variational calculations. Jauregui and Bunge [71] 
used a configuration-interaction(CI)-expansion of 108 configuration terms and obtained 
the energy -0.125 354 716 6 a.u.. They analyzed the convergence pattern of their 
computations and extrapolated to -0.125 355 08(10) a.u.. The largest Hylleraas-type 
computation for this state, done by Drake and then repeated by [73] gave only -0.125 
335 6 a.u., though they used about 50 thousand terms containing powers of r\ up to 82. 
The best result, as we know, was obtained by Bylicki and Bednarz [72] who applied the 
Hyllerass configuration-interaction correlated expansions and used 1442 terms to give 
-0.125 355 451 24 a.u.. More efforts are needed than that for the ground state. In our 
approach, this problem is reconsidered from a new viewpoint, and the finite mass of 
proton is naturally taken into account. 
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Figure 5. Radial density probability distribution P(r) for °°H~( 3 P), H(2p), H(3p) 
and H(4p). 

We compare our results with other calculations in table El The convergence is just 
as good in this case as for the ground state, and the leading term of all the sequences 
is estimated to be x 5 . Our energy eigenvalue is 0.125 355 451 242 a.u. for °°H~, better 
than most of previous calculations. 

As a negative ion, the extra electron is significantly affected by the other electron. 
FigureEJshows the distribution of n, r-ii ^3 and cos# of the °°H~. The radial distribution 
has a very long tail; even at r = 60 the probability does not vanish. If the electron- 
electron interaction is ignored so as to approximate this state by two 2p-state electrons, 
the charge distribution of the 2p state would to vanish at about r = 15 (see figure EJ, 
very far from r = 60. This long tail exponentially decays with a decay length I ~ 18.3, 
with a form r a exp (— r/l) with a ps 0.5. In this case, the independent-electron picture is 
very far removed from the present excited state. The variational method is observed to 
be less accurate for this state due to the very strong electron-electron correlation. For 
some of the expectation values, such as < r; >, comparable results can not be found 
in the literature. To our knowledge, our calculated result on this state is probably the 
most accurate so far. 

It is interesting to observe that in this state , the two electrons form nearly an angle 
of 7r/2 relative to the nucleus (the value of < cos(ri 2 ,ri 3 ) >^ 0). Also, the value of 
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< ri > is nearly twice that of the ground state. 
5.3. Helium and helium-like ions 

Except for the hydrogen atom, helium atom (or helium-like ions) is perhaps the simplest 
system in quantum mechanics. However, the calculation of its properties is not trivial 
and represents a real challenge. It has been 80 years since the beginning of quantum 
mechanics, and there are still efforts trying to understand the system from some new 
perspectives. 

In atomic theory, the shell model is the starting point to understand complex atoms. 
This model is an independent-electron picture, and the interaction between electrons 
are averaged by a mean field. For helium, we denote one of the electrons as the inner- 
electron, and the other as the outer-electron. But when the correlation of two electrons 
is strong, this picture fails. For example, in the shell model language the two electrons 
in the ground state are in the same orbit. Interaction would thus become important, 
which would make the independent-electron picture not suitable. For the doubly excited 
states in which both electrons are not in the ground state, interaction can affect the 
slow-moving electrons very significantly. 

5.3.1. The S states of helium These states have been extensively studied, and there 
are many methods used to treat the S states. For extrapolation, we estimate the leading 
term of E(x) (for all of the S states) to be x 3 , and a list of the limiting values are given 
in table El after performing the Richardson extrapolation. Except for the ground state 
1 S e , the other states are close to the independent-electron picture. In tables El and El 
some of expectation values are given for the 1 ' 3 S e states of °°He. 

5.3.2. The odd parity P states The leading term of E(x) of all of these states is 
estimated to be x 2 . Table |U shows the convergence for the three low-lying 1 > 3 P° states 
of 00 He. 

5.3.3. The parity-unfavorable 1 > 3 P e ; 1,3 _D° ; 1 > 3 F e states The doubly excited states 1 > 3 P e 
and 1 ' 3 D° have been observed in experiments long ago. Dolye et al j2Zj calculated the 
positions of the l ' z P e and 1,3 D° states of the helium isoelectronic sequence using the 1/Z 
expansion method. By using a variational scheme Bhatia [H] had calculated the 1,3 D° 
states of helium again and obtained a few higher-precision results. The best results for 
the states of helium were obtained by Goodson et al [TS], who took the advantage of 
the interdimensional degeneracies of the problem and used the variational method in 
the five-dimensional space to obtain the results. The observation of the 1,3 _F e states is 
rarely reported, and its energies had been computed by Galan and Bunge jZH]- In this 
subsection, we report the energies re-calculated by our approach. 

For the 1 > 3 P e states, we compare our results with other calculations in table El As a 
non- variational method, the performance of our approach is excellent. The distributions 




r r 



Figure 6. Distributions of r\, T2 (or r-&) of the five low-lying 3 P e states of °°He (ground 
state plus four excited states). The solid line is the distribution of T2 or r% (same), 
and the dashed line is the distribution of r\. For comparison, we also plot the radial 
density distribution of the neutral He(2p), H(2p) and H(3p) states in (5) (normalized 
to 1/2). 
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(X, Y) 




Figure 7. The coordinate system of • 



of ri, T2, r3 and COS6 1 are shown in figure El for °°He. Some properties of these states 
are given in tables [7| and |H1 It should be mentioned that only in the case of Ei of the 
1 P e state is our result higher in energy than the variational approach, but has more 
significant figures. We think that in this case our result is correct, because our excited 
states' eigen-energies and their related wavefunctions were all obtained simultaneously 
with the ground state properties. Hence the accuracy of our excited states' properties 
are insured by the demonstrated accuracy of the ground state eigenvalue. For helium-like 
ions {Z — 3 ~ 6), the eigen-energies are summarized in table fTUl 

We have also calculated four low-lying 1 ^D° states of °°He. We compare our results 
with Bathia's in table ^2 The three low-lying states 1 F e states of Li + ion are given in 
table H21 It is seen that our results are better than those obtained by other methods. 

5-4- H\ and Ps~ 

In order to obtain a better understanding of the underlying dynamical property of the 
Coulomb three-body system with two identical particles, three typical systems are the 
most interesting: (1) where the mass ratio m\ is large compared to m 2 , m 3 in the 
identical particles pair, e.g., for a helium atom or a helium-like ion (such as the H 
ion), (2) where mi = m 2 = m 3 for a positronium negative ion Ps~, and (3) where mi is 
very small, e.g., for a hydrogen molecular ion H^. It should be noted that the case of 
hydrogen molecular ion is very close to a one-body system in which the single electron 
moves in the field of two heavy protons. 

5.4-1- The hydrogen molecular ion H\ We have calculated 15 low-lying 5* states for 
this ion, the most achieved so far, and the numerical results are summarized in table!13l 
For exploring the properties of H^" , the coordinate system (X, Y, R) as shown in figure [7| 
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Figure 9. The function \^(X,Y,R = Rq)\ for the ground state of H^~, where 
Rq = 0.5, 1.0. The electron density is noted to peak directly above the positions 
of the protons, with a sharp cusp. 

is used, and the wave function has the form: 

V = V(X,Y,R). (152) 

Here R is the distance between the two protons, and X, Y denote the coordinate of 
the electron relative to the two protons. Figure |H1 shows the distribution of R for the 
calculated states. It is obvious that there is vibrational motion between the two protons, 
similar to an anharmonic oscillator. In this case the attraction between the two protons 
is clearly mediated by the electron. Figures l9l and fTOl show the function \"ff(X,Y, R)\ 
of the ground state when R = 0.5, 1.0, 1.5, 2.0. It is clear from these pictures that 
the electron probability density is the highest directly above each of the protons, with 
a sharp cusp at the proton locations. 

5.4-2. The positronium negative ion Ps~ The experimental and theoretical studies of 
the positronium negative ion Ps~, consisting of two electrons and one positron, have 
attracted considerable interest since the work of Wheeler [H3], who proved Ps~ to be 
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0.15 




Figure 10. Same as in figure E| but Rq = 1.5, 2.0. 

stable against dissociation into a free electron and a positronium atom. The formation 
of Ps~ has been discussed in the e + — He jHS] and Ps — H (HH] scattering calculations. The 
binding energy has been calculated variationally by several authors. Recent interests 
include the calculations of autoionizing doubly-excited states and the investigations 
of the possible existence of the so-called second bound state 3 P e . These positronium 
negative ions have been observed in the laboratory by Mills [HZ1IEE]- A large number of 
the doubly-excited states of a positronium negative ion were calculated by the method 
of complex-coordinate rotations |89| . 

We compare our results with other calculations in table El Figure ^2 shows the 
distributions for n, r 2 (or rs) and cos# between the two electrons. It is seen that the 
two electrons plus the positron form a triangle in which the angle formed by r 2 and 
r% is greater than 60 degrees, consistent with the fact that the separation between the 
two electrons is larger than their respective distances to the positron. It should be 
noted that the curve of the cos 6 distribution is very different from that of the °°H _ (see 
figure EJ). Since the only difference between the two systems lies in the mass of the 
positive charge, the comparison is useful to delineate the importance of mass ratio. In 




r 



Figure 11. Distributions of n, r2 (or r 3 ) and cos£? of the ground state of Ps~. The 
solid line is the distribution of r 2 or r 3 (same), and the dashed line is the distribution 
of fx. Inset is the distribution of cos 9. 

the present case, the mass of the positive charge is the same as that of the electron, so 
the positron is envisioned to be much more mobile, whereas in the case of H~ the proton 
is more localized. This comparison is perhaps helpful in answering the question [§0]: Is 
the positronium system a molecule or an atom? 

6. Further developments 

The kinetic energy operator approach is noted to have applications potential for a 
number of interesting problems. The approach is obviously applicable to the three- 
body problem in two dimensions. The addition of a perpendicular magnetic field plus 
a circularly symmetric potential is noted to be possible in this case, without breaking 
the symmetry that is essential to the solution approach. Thus it would be particularly 
interesting to examine the interaction of three electrons under a strong magnetic field, 
a configuration which might be relevant to a fractional quantum Hall state. 

The fact that we can reduce the dimensionality of the three-body problem also 
promises numerical efficiency in the case of three-body scattering, which requires very 
extensive computational resources at present [5]. There is also the possibility of carrying 
out time- dependent calculations, with implications to the enumeration of doubly excited 
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states in the continuum. 



Appendix A. Eigenfunctions of A® 

Appendix A. 1. Eigenfunctions of A® 

In this appendix, we give the details for computing the eigenfunctions of the operator 

+ 



i - r — , 

5 ^£2 £ Q£ < £2^2' 

and list some useful properties of the Jacobi polynomials. 

Let x = 2£ 2 - 1 and ip = e ±im/3 £ m F(x), where m = 0, 1, 2, 
computation shows that 

^(Jty = e ±im/3^m - m{l + m + 2)F(x)} , 

where 



x) 



Al(l-x 2 )F (x) 



m 



2 + - + m j x 



F(x) 



(A.l) 

Straightforward 
(A.2) 

(A.3) 



The Jacobi polynomial P^' b \x) satisfies the following differential equation \D := ^j: 

(1 - x 2 )D 2 Pt b \x) + [b-a-(2 + a + b)x] DP^ b \x) 

+n(l + a + b + n)P^ b) (x) = 0. (A.4) 

The orthogonal relations of Jacobi polynomials |P^ a,fe ^(x), n = 0, 1, 2, • • -j are given by 

/ (1 _ ir(1 + ^.(^wd, = f:r^ + a 2 vt ( : + :::' ^.(a.5) 

7-i n!(l + a + o + 2n)i (1 + a + b + n) 

From ()A.3|) and ()A.4|) . we can see that the following set of functions constitute a complete 
family of eigenfunctions of A( l \ i.e., 



41,n(£,/3) = e ±im ^ m Pn^ l \2e - 1), m,n > 



with their eigenvalues given by 



A (0 



m,n 



4n(l + - + m + nj + m(7 + m + 2) 



(A.6) 



(A.7) 



The orthogonal relations become 

-1 4U41n 2 £(l-^ /2 d£d/3 



7T JO 



7rr(l + | + n)r(l + m + n) 
ra!(l + | + m + 2n)r(l + \ + m + n) 



(A. 
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Appendix A. 2. Recurrence relations of the Jacobi polynomials and the matrices of B\ 
and B 2 

Appendix A. 2.1. Matrices of B\ and B 2 For the two linear differential operators 

_ „d sm(3 d ,7-, ■ n 9 cos (3 d 

i3i=cosp— — — j an d B 2 = smp — + 



we have 
(i) The case of m = 0: 
r (0 



S!J l i = (2ie^ + 2ie-^)M 



.7' 



£ 2 jW = (-2ie^ + 2ie-^)^" 0) («)- 
(ii) The case of positive indices: 



+ e 



i(m— l)/3^m— 1 



mP!' ,m) (i) + (1 + x)DP& m \x) 



BoJ i £' r 



2ie i{m+1 ^C +1 DPn"' m \x) 



+ ie 



i(m— l)/3^m— 1 



mP^' m) (i) + (l + x)rf' m) (i) 



(in) The case of negative indices: 



(0 



B x J 



R T {1) 

n 2'J-m,n 



2e 



_|_ g-i(m-l)j9twi— 1 



m# m) W + (l + ^ ,m) (i) 



2ie- i (^+i)^^+i J DP^' m) ( x ) 



— le 



-i(m— l)/3^-m— 1 



mP^' m) (x) + (1 + i)^' m) (i) 



Therefore, we need explicit formulas which express DP> m) (x) (mP> m) (x) + (1 + 
x)DPn 2 ' m \x)) as linear combination of Pl 2 ' m+1 \x)(P^ 2 ' m 1 (x)). Here, the existence 
of a collection of recurrence relations among Jacobi polynomials provides a handy tool 
for such a purpose. 



Appendix A. 2. 2. Recurrence relations of the Jacobi polynomials We recall some of the 
mixed-type ones from §138 of Rainville |94j : 



DP^ b \ 



x 



(l + x)DPi a 



x 



-{l + a + b + n)Pt\ l ' b+1 \x) 

1 -(b + n)Pt\^ b \x) + i( fl + n)P^ +l \x) 
nPi a ' b \x) + (b + n)P^ h \x) 



(a + b + 2n)P^ h - 1 \x) =(a + b + n)P^ b \x) + (a + n)P { ^(x) 
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(a + b + 2n)P^ 1 ' b) (x) =(a + b + n)P^ b) (x) - (b + n)P^ (x) 
P^\x)-P^ b \x) = P^(x). 

We have the following specific recurrence relations: 



DP^ b) \ 



x) 



-(a + b) + n 



P^\ +1 \x) 



b + n 
a + b + n 



DP£${x 



bP^ b \x) + (1 + x)DP^ b \x)] = {b + n)P^ a+1 ' b - x \x) 



a + o + n 1 



+ 



[X . 



a + 6 + n 

By using the above formulas, we get the explicit expressions for the matrices of B\ and 
B 2 . The results are given in section 3. 



Appendix B. IMT Integration Scheme 

Consider the real value function f(x) defined on the interval [0, 1], which is continuous 
and differentiable sufficiently many times on [0,1]. For the N-point trapezoidal rule, 
there is the Euler-Maclaurin formula 



£/(*)dx=i 



N-l 



^(/(0) + /(l))+£/(n/AO 



m— 1 d i 

(-i;- 

r=l 



+ E 



/(2 ,-l) (1) _ /( ,r.-I, (0) + ^ (R1 



?(2r-l) 



(2r)!iV 2r 

where _B r is the r-th Bernoulli number and R m is the reminder term. Noting that the 
expression for the truncation error depends only on values of the functional derivatives 
at the integration end points, we design a change of variable, so that 

/(2r-i)(l) = /(2r-i)(0) for r = 1, 2, • • ■ . (B.2) 

It follows from (IB~T| that all error terms will vanish, and high precision can be achieved. 
The IMT- rule j^SJ ESI) a ^ so denoted the "double-exponential formula", is one type of 
integration technique. Consider 

1 



Q 



1 ( 1 

exp — 



dt 



/o v t l-t, 
0.00702 98584 06609 65623 92412 70530 



Define the function ip(t) by 
<p(t) 



t / i 

ip(r)dr, (f'(t) = —exp 
Q 



(B.3) 



(B.4) 



t l-t, 

The transformation x = ip(t) maps the variable x G [0, 1] into the variable t G [0,1], 
and we have a new expression for the integration of f(x) over [0,1]: 



f(x)dx = g(t)dt, with g{t) = f{<p(t))<p'(t). 



(B.5) 
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If the function f(x) is differentiable infinitely many times on (0, 1) and has an algebraic 
singularity such as x a or (1 — Xf > —1), then the function g(t) also can be 

differentiable many times and all its derivatives vanish at the two ends of the interval 
[0, 1] due to the strong singularity of <p'(t), so that 

^M(O) = g {m) (l) = 0, for m = 0,1,2,- ••. (B.6) 

Applying the trapezoidal rule on (|B.5|) . we get the approximation for the integration: 

i N-l 



s 



N 



N 



£ wWfixW), and X W = tp(n/N),wW = tf/(n/N). 



n=l 



If we define 



then obviously 



and 



9(t) 



g(t) exp(i27rfci) dt, 

oo 

y~] Cfc exp(— i2nkt) 



k=—oo 



I f{x)dx = [ 1 g(t)dt = c . 
Jo Jo 



The approximation Sn is given by 

N-l 

Sn 



Y JV— l oo 

— 9(n/N) = c + Y.( c pN 

iV n=0 p=l 



-pN 



)• 



The truncation error En is estinated as: 



En — S 



N 







git) dt 



— Yj( C pN + C- p n) 
p=l 

oo 

= 2^Re c pN 
p=i 

~ 2Re c^. 

For f(x) with singularity x a or (1 — x) a , an estimate of e(N, a) is given by 
v 7 ^ (a + l)V4+a 



(eQ) a+1 
• exp 



(2vriV) 3 /^ 



^/47r(a + 1)N ■ cos ^/47r(a + ljiV + 



3 + 4a 
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(B.7) 
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(B.9) 



(B.10) 
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Table 1. Comparison of the ground state of °°H and H (m p = 1836.152701) with 
other theoretical calculations. 
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503 


b 





.310 


815 


007 


66 c 


< ?*2 > = < T3 > 


2 


,710 


178 


27 




2 


.712 


095 


6 






2 


,710 


178 


278 


34 a 


2 


.712 


095 


626 


51 a 




2 


,710 


178 


263 


b 


2 


.712 


095 


621 


4 C 


< n > 


4. 


,412 


694 


50 




4 


.415 


692 


6 






4 


,412 


694 


497 


79 a 


4 


.415 


692 


603 


31 a 














4 


.415 


692 


593 


4 C 


< r\ >=< r\ > 


11 


.913 


699 


6 




11 


.931 


747 


7 






11 


,913 


699 


681 


6 a 


11 


.931 


747 


760 


a 




11 


,913 


699 


235 


b 


11 


.931 


747 


62 c 




<r\> 


25 


.202 


025 


2 




25 


.237 


175 








25 


.202 


025 


298 


2 a 


25 


.237 


175 


614 


i a 














25 


.237 


175 


34 c 




< cos(r 12 ,r 13 ) > 


-0. 


,105 


147 


693 


7 


-0 


.104 


996 


606 






-0 


,105 


147 


693 


566 a 


-0 


.104 


996 


606 


303 l a 


< cos(r 12 ,r 32 ) > 





649 


871 


581 







,694 


795 


647 









649 


871 


581 


193 9 a 





.694 


795 


646 


586 a 



a Rcf. HD] b Ref. [nil c Ref. ^ 
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Table 2. Comparison of the 3 P e state of °°H~ and H~ (m p = 1836.152701) with 
other theoretical calculations. 







H 


LP 
-III 


n ioc; qc.c /tc.i O/io 
U.lZO oOO 401 Z4Z 


H IOC OfiQ 1 C.7 HQ/1 
U.IZO Zoo 10/ Uo4 




n i oc; o C/l 7a 
U.lZO o04 / 






U.IZO oOO 401 Z4 






n 1 oc. qc./i 7nc. c 

U.IZO o04 /UO 


n 1 oc. oqo ^01 c 
u.izo zoz oyi y 




A IOC qci od 

U.IZO OOl 


n 10c 070 n d 

u.izo z < y u 




fl IOC Q C C HQ G 

U.IZO 000 Uo 




< t/^2 > — < 1/^3 > 


n 1 fin con S7fi /i 

U. 1DU OZU / 4 


n 1 fin ^07 1 fi c. 
u.iou oy 1 loo 


< 1/n > 


0.0703 308 548 


0.0702 280 17 




0.0706 30 a 




< r 2 >=< r 3 > 


11.657 657 7 


11.683 161 8 


< ri > 


19.585 091 


19.632 096 8 




19.237 a 




< >=< r§ > 


271.263 4 


273.015 9 


< r\ > 


557.259 2 


560.744 




517.09 a 




< cos(ri2,ri 3 ) > 


-0.093 867 209 4 


-0.093 625 370 


< cos(ri2,r 32 ) > 


0.651 280 257 


0.651 086 687 


a Ref. [73, b Ref. [721, CRef - EE dRef - ES1> eRef - EH 



Table 3. Comparison for the energies of the five low-lying 1,3 5 e states of °°He and 
He with others calculations. 





-E^S") -E( 3 S e ) 


This work 


Ref. QUI This work 


Ref. QUI 




00 He 




2.903 724 377 03 


2.903 724 377 034 119 2.175 229 378 2 


2.175 229 378 236 


2.145 974 046 


2.145 974 046 054 2.068 689 067 


2.068 689 067 47 


2.061 271 98 


2.061 271 989 7 2.036 512 


2.036 512 083 


2.033 586 


2.033 586 7 2.022 6 


2.022 618 


2.021 17 


2.021 17 2.015 


2.015 377 




He (m a = 7294.299507) 




2.903 304 557 7 


2.174 930 190 6 




2.145 678 587 


2.068 405 243 




2.060 989 07 


2.036 232 73 




2.033 307 6 


2.022 280 




2.020 77 


2.013 34 
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Table 4. Convergence study for the three low-lying states of °°He. 



(M,N,P) -E -E x -E 2 















1 po 


















(290,160,60) 


2. 


.121 


781 


471 


934 


2.045 


206 


888 


028 


2 


,007 


282 


642 


688 


(310,170,60) 


2. 


.122 


001 


335 


701 


2.046 


159 


516 


026 


2 


.009 


274 


234 


036 


(330,180,60) 


2. 


.122 


187 


837 


151 


2.046 


981 


326 


187 


2 


Oil 


018 


748 


522 


(350,190,60) 


2. 


.122 


347 


396 


423 


2.047 


695 


333 


577 


2 


,012 


556 


633 


845 


(370,200,60) 


2. 


.122 


484 


960 


882 


2.048 


319 


671 


610 


2 


,013 


920 


175 


194 


(390,210,60) 


2. 


.122 


604 


392 


689 


2.048 


868 


793 


171 


2 


,015 


135 


435 


817 


(410,220,60) 


2. 


.122 


708 


742 


385 


2.049 


354 


337 


962 


2 


016 


223 


670 


874 


(430,230,60) 


2. 


.122 


800 


445 


240 


2.049 


785 


768 


826 


2 


,017 


202 


374 


649 


Extrap. 


2. 


123 


842 


8 




2.055 


149 






2 


,031 









Ref. [SOI 


2. 


123 


843 


085 


800 


2.055 

3 po 


146 


355 


4 


2 


,031 


069 


591 




(290,160,60) 


2. 


131 


433 


366 


926 


2.049 


022 


358 


625 


2 


,009 


815 


747 


062 


(310,170,60) 


2 


131 


618 


811 


233 


2.049 


898 


916 


341 


2 


Oil 


718 


260 


682 


(330,180,60) 


2. 


131 


775 


987 


950 


2.050 


654 


119 


626 


2 


,013 


383 


329 


203 


(350,190,60) 


2. 


131 


910 


360 


569 


2.051 


309 


458 


014 


2 


,014 


849 


959 


935 


(370,200,60) 


2. 


.132 


026 


134 


097 


2.051 


881 


840 


210 


2 


,016 


149 


272 


454 


(390,210,60) 


2. 


.132 


126 


587 


511 


2.052 


384 


725 


339 


2 


,017 


306 


376 


122 


(410,220,60) 


2. 


132 


214 


308 


072 


2.052 


828 


939 


042 


2 


,018 


341 


737 


893 


(430,230,60) 


2. 


132 


291 


359 


280 


2.053 


223 


271 


537 


2 


,019 


272 


194 


566 


Extrap. 


2. 


133 


164 


06 




2.058 


081 






2 


,032 


36 






Ref. [HOI 


2. 


133 


164 


190 


534 


2.058 


081 


081 


6 


2 


,032 


324 


325 





Table 5. Properties of the five low-lying 1 S e states of °°He. 



State 


1 l S e 


2 1 S e 


3 1 S e 


4 1 S e 


5 1 S e 


< l/r 2 > 


1.688 316 800 71 


1.135 407 686 


1.058 514 75 


1.032 484 8 


1.021 298 


< 1/ri > 


0.945 818 448 80 


0.249 682 652 


0.111 514 95 


0.062 760 2 


0.041 429 


<r 2 > 


0.929 472 294 87 


2.973 061 12 


6.511 676 


11.550 6 


17.573 


<n> 


1.422 070 255 5 


5.269 696 20 


12.304 521 


22.368 1 


34.407 


<r| > 


1.193 482 995 


16.089 233 2 


85.890 18 


281.33 


660.7 


<r{> 


2.516 439 312 9 


32.302 380 3 


171.838 66 


562.71 


132 1.4 


< cos(ri2,ri 3 ) > 


-0.064 202 614 217 


-0.014 657 043 3 


-0.004 317 036 7 


-0.001 795 63 


-0.000 958 99 


< cos(ri2,r 32 ) > 


0.648 017 667 47 


0.557 144 578 


0.526 466 53 


0.515 120 


0.510 057 
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Table 6. Properties of the five low-lying 3 S e states of °°He. 



State 


2 3 S e 


3 3 S e 


4 3 S e 


5 3 S e 


6 3 S e 


< l/r 2 > 


1.154 663 198 4 


1.063 661 050 


1.034 491 5 


1.021 740 


1.019 08 


< i/n > 


0.268 197 633 6 


0.117 316 73 


0.065 253 6 


0.042 086 


0.037 343 


< r 2 > 


2.550 464 78 


5.856 041 6 


10.661 70 


16.702 8 


19.945 


< n > 


4.447 538 89 


10.998 910 1 


20.592 6 


32.666 8 


39.150 


< r\ > 


11.464 340 5 


68.710 08 


238.596 


596.06 


870.56 


< r\ > 


23.046 235 5 


137.478 44 


477.225 


119 2.20 


174 1.3 


< cos(ri2,ri 3 ) > 


-0.015 839 217 1 


-0.004 245 085 8 


-0.001 686 96 


-0.000 853 68 


0.000 753 55 


< cos(ri2,r 32 ) > 


0.562 788 947 4 


0.528 301 29 


0.515 917 


0.510 32 


0.509 16 



Table 7. Properties of the five low-lying 1 P e states of °°He. 



State 


3 1 P e 


4 i pe 


5 1 P e 


Qlpe 


7 lpe 


< l/r 2 > 


.320101054072 


.28662333040 


.27259779205 


.2653459435 


.261119621 


< > 


.119911271506 


.066410151506 


.0455874718 


.0320333447 


.021194314 


< r 2 > 


5.68748370235 


10.289515048 


16.39163259 


23.9930126 


33.043005 


< rt > 


9.3831001400 


18.33321725 


30.43563233 


45.5874716 


63.658 480 


<r\> 


48.197972382 


187.5444301 


517.8278870 


1160.99752 


2260.9388 


<r\> 


97.935626136 


376.133414 


1036.36343 


2322.49992 


4522.255 


< cos(ri2,ri 3 ) > 


-0.031985135411 


-0.01274515627 


-0.00621248298 


-0.0034696118 


-0.002133103 


< cos(ri2,r 32 ) > 


.608746828961 


.56272353076 


.54033973838 


.5280463517 


.520635797 



Table 8. Properties of the five low-lying 3 P e states of °°He. 



State 


2 3 P e 


3 3 P e 


4 3pe 


5 3 P e 


6 3 P e 


< l/r 2 > 


.41809816667770 


.312484187824 


.28403411128 


.27140054988 


.2646919161 


< 1/ri > 


.2513923553540 


.114310953848 


.06440206759 


.04109304812 


.028447092 


< r 2 > 


3.089879033146 


6.460157023 


11.30744532 


17.67087220 


25.5392910 


< n > 


4.676371886388 


10.831719973 


20.33743461 


32.98000076 


48.6725310 


<r\> 


11.79098787761 


65.49966426 


232.2699200 


609.9032260 


1325.7638 


< r\ > 


25.06816973162 


132.73779825 


465.6256393 


1220.528804 


2652.0392 


< cos(ri2,ri 3 ) > 


-0.0714037175082 


-0.031229607240 


-0.01339606321 


-0.00680633501 


-0.0038988799 


< cos(ri2,r 32 ) > 


.6800429403893 


.603424804436 


.56048920421 


.53923906130 


.527421792 
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Table 9. Comparison for the energies of the five low-lying 1 ' 3 P e states of °°He and 
He with other calculations. 





-E^P") 




-E( 3 P e ) 


This work 


other results 


This work 


other results 


00 He 


0.580 246 472 594 


0.580 246 472 594 392 a 


0.710 500 155 678 3 


0.710 500 155 678 334 3 a 




0.580 246 472 594 388 b 




0.710 500 155 678 23 b 








0.710 500 155 656 78 c 


540 041 590 93 


540 041 590 938 1 a 


567 81 2 898 725 1 


567 812 898 725 31 a 




0.540 041 590 938 52 b 




0.567 812 898 725 16 b 


0.524 178 981 8 


0.524 179 01 a 


0.535 867 188 767 


0.535 867 188 71 a 


0.516 208 609 4 


0.516 03 a 


0.522 254 575 706 


0.522 253 a 


0.511 624 834 




0.515 160 203 83 






He (m a 


= 7294.299507) 




0.580 165 768 725 


0.580 165 768 308 4 b 


0.710 396 457 557 


0.710 396 457 021 81 b 




0.580 165 768 d 




0.710 396 457 d 


0.539 967 178 01 


0.539 967 177 633 b 


0.567 733 870 122 


0.567 733 869 714 03 b 




0.539 967 2 d 




0.567 733 87 d 


0.524 106 954 1 




0.535 793 284 74 




0.516 137 755 4 




0.522 182 770 7 




0.511 554 646 4 




0.515 089 462 





a Ref. EH, b Ref. EU, c Ref. [SU, d Ref. [S2]. 



Table 10. The eigen- energies for the 1 > 3 P e states of helium-like ions (Z — 3 ~ 6). 







-£;( 1 P e ) 






-P( 3 P e ) 








- J B( 1 P e ) 








P( 3 P e ) 












Z 


=3 
















Z 


=4 










1 


.401 


410 


927 


020 


1.796 


648 


099 


720 


2. 


583 


994 


187 


432 


3 


.382 


712 


420 


777 


1 


.269 


787 


972 


287 


1.373 


589 


535 


176 


2. 


312 


232 


549 


880 


2 


.540 


768 


798 


391 


1 


.214 


520 


958 


34 


1.260 


545 


265 


63 


2. 


.194 


982 


661 


13 


2 


.298 


106 


658 


34 


1 


.185 


881 


076 


3 


1.210 


287 


247 


90 


2. 


133 


407 


126 





2 


.188 


554 


817 


30 


1 


.169 


099 


559 




1.183 


580 


072 


8 


2. 


.097 


031 


380 




2 


.129 


923 


592 


6 










Z: 


=5 
















Z 


=6 










4 


.127 


776 


355 


386 


5.468 


730 


984 


923 


6, 


032 


706 


408 


00 


8 


.054 


724 


273 


276 


3 


.667 


237 


481 


32 


4.069 


185 


038 


025 


5, 


.334 


768 


645 


76 


5 


.958 


774 


148 


469 


3 


.465 


482 


470 


00 


3.648 


310 


490 


768 


5, 


.025 


999 


284 


37 


5 


.311 


076 


536 


74 


3 


.358 


735 


604 


6 


3.456 


902 


285 


67 


4. 


861 


852 


996 


5 


5 


.015 


282 


040 


9 


3 


.295 


388 


021 


8 


3.354 


092 


319 





4. 


.764 


160 


404 




4 


.856 


057 


141 
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Table 11. Comparison for the energies of the five low-lying 1 ' 3 £)° states of °°He with 
Bhatia's [711 



-EfE 


,0) 


-E( 3 D°) 


This work 


Bhatia 


This work 


Bhatia 


0.563 800 420 4 
0.534 576 384 
0.521 659 00 
0.514 833 4 


0.563 800 405 
0.534 576 015 
0.521 642 77 
0.514 269 06 


0.559 328 263 
0.532 678 600 
0.520 703 44 
0.514 288 2 


0.559 328 25 
0.532 678 075 
0.520 693 865 
0.514 235 78 



Table 12. Convergence study for the three low-lying 1 F e states of Li + ion. 



(M,N,P) 


-E a 




-E l 




(310,170,40) 


1.252 443 716 475 54 


1 


.206 193 104 449 


1.181 010 621 655 


(330,180,40) 


1.252 445 059 532 88 


1 


.206 204 195 186 


1.181 059 878 109 


(350,190,40) 


1.252 446 091 667 05 


1 


.206 212 775 544 


1.181 098 378 639 


(370,200,40) 


1.252 446 895 466 84 


1 


.206 219 497 607 


1.181 128 811 337 


(390,210,40) 


1.252 447 529 007 27 


1 


.206 224 824 200 


1.181 153 115 955 


Extrap. 


1.252 450 636 


1 


.206 251 57 


1.181 279 1 




1.252 445 1 Ref. 










1.252 258 Ref. [SS| 









Table 13. 15 low-lying S states of 



11 


—E n 


n 


—E n 


n 


-E n 





0.597 139 


5 


0.552 841 


10 


0.521 699 


1 


0.587 156 


6 


0.545 593 


11 


0.517 002 


2 


0.577 752 


7 


0.538 858 


12 


0.512 827 


3 


0.568 909 


8 


0.532 631 


13 


0.509 189 


4 


0.560 609 


9 


0.526 911 


14 


0.506 11 
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Table 14. Comparison of the ground state of Ps with other theoretical calculations. 



-E 


< l/r 2 >=< l/r-3 > 


< 1/ri > 


0.262 005 070 2 


0.339 821 023 


0.155 631 905 7 


0.262 004 857 a 


0.339 831 3 a 


0.155 654 3 a 


0.262 005 070 b 


0.339 821 02 b 


0.155 631 90 b 


0.262 005 070 232 94 c 


0.339 821 023 06 c 


0.155 631 905 653 c 


0.262 005 070 232 978 d 


0.339 821 023 059 27 d 


0.155 631 905 652 66 d 


< r 2 >=< r 3 > 


< n > 


< ^ > = < J3 > 


5.489 633 2 


8.548 580 6 


48.418 936 


5.488 352 a 


8.546 111 29 a 


48.379 317 a 


5.489 633 3 b 


8.548 580 8 b 


48.418 936 b 


5.489 633 252 c 


8.548 580 655 c 


48.418 937 2 C 


5.489 633 252 38 d 


8.548 580 655 16 d 


48.418 937 230 d 


< r\ > 


< cos(ri2,ri 3 ) > 


< cos(ri2,r 32 ) > 


93.178 63 


0.019 769 632 8 


0.591 981 70 


93.100 697 a 






93.100 633 b 






93.178 633 80 c 






93.178 633 855 d 


0.019 769 632 816 7 d 


0.591 981 701 149 2 d 



Ref. b Ref. [HH c Ref. ^ d Ref. ^ 



